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ABSTRACT 


Inversion  techniques  are  used  to  infer  bottom 
sedimentary  characteristics  using  AN/SQS-53C  Reverberation 
Level  (RL)  data  gathered  during  Exercise  LWAD  99-1. 
Analysis  was  conducted  to  determine  the  correlation  between 
sediment  type  and  beam  RL  time  series  statistics,  including 
slope  and  detrended  higher  order  statistics  (standard 
deviation,  skew,  and  kurtosis) .  A  geographic  contour  plot 
of  the  spatial  distribution  of  sediments  was  generated  from 
the  deviation  of  the  RL  slope  of  individual  beamform  data 
from  the  area-wide  average  RL  slope.  The  resulting  plot 
shows  features  similar  to  that  of  the  pre-exercise 
sedimentary  characterization  map.  Higher  order  statistics 
are  found  to  be  non-Gaussian  which  has  signal  processing 
implications . 
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I .    INTRODUCTION 

A.    OBJECTIVE 

The  objective  of  this  thesis  is  to  develop  a  real  time, 
high-resolution  method  for  determining  bottom  sediment 
acoustic  properties  in  a  shallow  water  environment.  This 
study  sought  to  determine  those  properties  through  the 
analysis  of  the  slope  characteristics  of  beam  reverberation 
data  and  higher  order  statistics  (HOS)  related  to  the  RL 
beam  time  series. 

Current  fleet  anti-submarine  warfare  (ASW)  operations 
are  limited  in  shallow  water  in  part  because  of  their 
dependence  on  bottom  sediment  type  maps  that  are  based  on 
large  area  averages  (low  resolution)  derived  from  few 
samples  [Ref .  1]  .  Such  maps  provide  limited  or  no 
information  on  the  bottom  roughness  or  the  sediment 
composition  beneath  the  water- sediment  interface.  They 
further  limit  the  accuracy  of  ASW  operations  because  of  the 
high  spatial  variability  of  the  sedimentary  composition, 
both  laterally  and  vertically,  that  is  a  common  feature  of 
shallow  water  continental  shelf  areas. 

This  thesis  examines  the  possibility  of  using  inversion 
techniques  to  deduce  sedimentary  properties  from  beam 
reverberation  level  (RL)  data.  The  use  of  RL  data  is 
practical  not  only  because  a  large  amount  of  data  is 


available  for  research,  but  because  the  RL  data  is  gathered 
by  active  fleet  assets  rather  than  an  oceanographic  research 
vessel.  In  this  study,  the  RL  data  was  gathered  during  the 
Littoral  Warfare  Advanced  Development  (LWAD)  Exercise  99-1, 
aboard  a  Ticonderoga  class  cruiser,  using  the  AN/SQS-53C 
tactical  sonar. 

B.    SHALLOW  WATER  ASW  PROBLEM 

Since  the  end  of  the  cold  war,  anti-submarine  warfare 
(ASW)  operations  have  moved  geographically  from  the  deep, 
open  ocean  environment,  to  the  near shore,  shallow  water, 
littoral  environment.  Cold  war  ASW  operations  had  their 
primary  emphasis  on  locating  and  tracking  Soviet  ballistic 
missile  submarines  (SSBN)  on  patrol  in  the  open  ocean.  Over 
the  last  10  years,  significant  geo-political  events  have 
occurred  which  have  shifted  ASW  emphasis  to  the  littoral 
environment.  With  the  downfall  of  the  Soviet  Union  came  a 
significant  decrease  in  their  SSBN  patrols,  while 
concurrently  new  threats  arose  in  such  places  as  the  Middle 
East  and  the  Korean  peninsula.  Smaller  countries  looked  to 
acquire  diesel-powered  submarines  for  coastal  defense. 
These  vessels  have  become  relatively  affordable  causing 
their  numbers  and  countries  of  ownership  to  proliferate. 
The  U.S.  Navy  ASW  doctrine  has  responded  to  this  paradigm 
change  by  placing  increased  emphasis  on  shallow  water 
operations  [Ref.  2]. 


Shallow  water  ASW  operations  are  significantly  more 
challenging  than  those  conducted  in  the  deep,  open  ocean. 
The  littoral  environment  includes  ocean  fronts  and  eddies, 
stronger  currents,  extensive,  often  steep  topographic 
changes,  increased  sea  life,  varied  bottom  sediment 
distributions,  and  increased  shipping  traffic.  These 
factors  combine  to  create  a  challenging  ASW  environment. 

C .    REVERBERATION 

One  of  the  most  significant  factors  contributing  to  the 
complexity  of  shallow  water  ASW  operations  is  the  spatial 
variability  of  the  sea  bottom  composition.  The  reflective 
nature  of  the  overlying  sediment  exerts  a  significant  effect 
on  active  sonar  reverberation  levels  (RL) .  Reverberation  is 
the  sum  of  all  random  scattering  mechanisms  in  the  ocean 
[Ref .  3,  p.  237]  .  In  deep  water  (Figure  1)  reverberation 
primarily  exists  as  a  combination  of  two  scattering 
mechanisms,  volume  scattering  and  sea  surface  scattering. 
Volume  scattering  occurs  due  to  scattering  from  marine  life 
within  the  water  column  and  from  inhomogeneities  in  the 
ocean  thermal /density  structure.  Sea  surface  scattering  is 
due  to  random  reflections  (scattering)  from  a  roughened  air- 
ocean  interface,  and  is  a  function  of  wind  speed  or  wave 
height . 

In  shallow  water  (Figure  2)  reverberation  from  the  sea 
bottom  typically  dominates  that  due  to  volume  or  sea  surface 


reverberation.  The  contribution  of  bottom  scattering  to  the 
total  reverberation  level  (RL)  is  often  so  severe  that  a 
target  echo  that  may  have  been  detected  in  deep  water  can 
not  be  detected  in  shallow  water  (Figures  1  and  2) . 

D.    BOTTOM  BACKSCATTERING 

Unlike  RL,  which  is  scattering  in  all  azimuthal 
directions,  backscattering  is  only  that  component  of  the 
scattered  energy  that  returns  to  the  receiver  in  monostatic 
(co-located  transmitter/receiver)  systems.  The  intensity  of 
the  backscattered  energy  from  the  sea  floor  that  arrives  at 
the   sonar   receiver   is   parameterized  as   the   scattering 

strength,  s  =10  loq  ^cat  [Ref.  3,  pp.  238-239].   Iscat  is  the 

b         I. 

inc 

intensity  of  the  scattered  sound  per  unit  area  at  a  unit 
distance  and  I.    is   the   incident  plane  wave   intensity 

inc  ■•-  -* 

impinging  on  the  scattering  area.  Measurements  of  RL  and 
Sb,  arising  from  sea  bed  backscattered  energy,  have 
indicated  that  Sb  is  a  function  of  grazing  angle  (9)  and 

bottom  type.  Backscattering  strength  traditionally  is 
represented  using  a  Lambert's  law  relationship  where  Sb  =  10 
log  |X  +  10  log  sin2(0)  .    0  is  the  grazing  angle  and  u 

represents  the  degree  of  reflectivity  of  the  sediment 
interface.  Based  on  the  deep  water  data  of  Mackenzie  [Ref. 
4]  and  others,  Urick  [Ref.  3,  pp.  271-273]  reports  that  over 


moderately  rough  topography  of  varying  sedimentary  makeup 
and  small  to  intermediate  grazing  angles  Mackenzie's 
constant   (k  =  10   log  (I)   averaged  about  -29  dB,   being 

essentially  invariant  in  frequency  and  grazing  angle.  More 
recently,  shallow  water  experiments  have  indicated  an 
association  of  Mackenzie's  constant  to  the  bottom  sediment 
type  and  that  the  grazing  angle  dependence  may  be  scaled  as 
sin(0)    over   muddy,    absorptive   bottoms   where   volume 

reverberation  from  the  upper  layers  of  the  sediment  may  be 
important  [Ref.  5].  This  finding  suggests  a  direct 
correlation  exists  between  backscattering  strength  and 
sediment  type.  With  this  knowledge,  it  is  possible  to  use 
inverse  techniques  to  infer  the  nature  of  the  sediment  type 
as  a  function  of  its  reflectivity  from  RL  data. 

E.    INVERSION  TECHNIQUES 

The  "forward  problem"  in  acoustics  are  those  methods  in 
which  one  predicts  the  characteristics  of  the  acoustic 
propagation  field  from  a  knowledge  of  the  governing  oceanic 
parameters  (e.g.,  sound  speed  structure,  sea  surface 
roughness,  geo-acoustic  properties,  etc.)  through  which  the 
sound  is  travelling.  "Inverse"  (or  inversion)  techniques 
are  the  exact  opposite  such  that  one  infers  one  or  more 
oceanic  parameters  through  the  analysis  of  acoustic 
propagation  data.  [Ref.  6] 


Inverse  techniques  have  a  wide  variety  of  applications 
including  high  resolution  geo-acoustic  modeling  [Ref.  7]  and 
simulation  of  broadband  signal  propagation  [Ref.  8]. 

Inversion  techniques  (IT)  are  used  in  this  study  to 
identify  features  in  the  beam  reverberation  time  series  from 
which  bottom  sedimentary  properties  can  be  deduced.  RL 
features  such  as  the  slope  of  the  RL  time  series  and  higher 
order  statistics  (HOS)  describing  the  shape  of  the  RL  time 
series  are  analyzed  to  identify  correlations  to  sediment 
type. 

F.    EXERCISE  LWAD  99-1 

The  Littoral  Warfare  Advanced  Development  (LWAD) 
Exercise  99-1  was  conducted  from  8  to  15  February  1999  in 
the  Gulf  of  Mexico  approximately  2  00  km  WNW  of  Key  West, 
Florida,  on  the  southwest  Florida  continental  shelf  (Figure 
3)  .  This  test  was  the  ninth  in  a  series  of  exercises 
performed  as  part  of  the  LWAD  project,  a  program  of  shallow 
water  exercises  designed  to  identify  new  technologies  which 
have  a  potential  to  be  transitioned  to  fleet  systems  in  the 
near  future.  The  LWAD  program  efforts  are  concentrated  on 
technologies  that  will  aid  in  the  location  and 
classification  of  diesel-electric  submarines  and  mines  in  a 
shallow  water  littoral  environment.  [Ref.  9] 

LWAD  99-1  included  several  experiments,  one  of  which 
was  for  the  continued  development  of  inversion  techniques. 


The  inversion  technique  experiment  had  the  following 
objectives:  (1)  record  AN/SQS-53C  data  in  a  geographically 
variable  environment  and  (2)  collect  ground  truth 
information  in  order  to  test  monostatic  and  bistatic 
inversion  techniques  [Ref.  10]. 

The  data  examined  in  this  thesis  was  collected  during 
event  31  conducted  on  9  February  1999  from  0730Z  to  1800Z  at 
LWAD  99-1  site  1  (Figure  3) .  The  USS  Vicksburg  ship's  track 
during  this  event  is  shown  in  Figure  4;  over  900  pings  from 
the  AN/SQS-53C  tactical •  sonar  were  recorded  during  the 
event.  Event  31  was  conducted  with  a  target  present;  the 
target  track  history  is  also  shown  in  Figure  4 . 

1.  Geo-Acoustics  of  the  South  Florida  Platform 
The  LWAD  99-1  exercise  area  lies  on  the  southwest 
portion  of  the  continental  shelf  bordering  the  west  coast  of 
Florida.  The  shelf  is  composed  of  a  karstic  Miocene 
platform  that  is  under  Pliocene-Pleistocene  and  Holocene 
sediments  [Ref.  11].  The  composition  of  the  surface 
sediment  was  analyzed  by  Holmes  from  a  variety  of  sources 
including  sediment  core  samples  and  seismic  data.  Figure  5 
is  a  depiction  of  the  surface  sedimentary  distribution  at 
LWAD  99-1  site  1  based  upon  Holmes'  analysis  and  represents 
the  most  definitive  characterization  of  the  upper  ocean 
sediment  composition  prior  to  the  commencement  of  LWAD  99-1. 


Howell  Hook  reef,  lying  along  the  western  portion  of 
site  1  at  a  depth  of  150  m  is  a  natural  border  between  the 
inner  shelf  and  the  shelf  slope.  Howell  Hook  is  composed  of 
coral  reef  material,  while  immediately  to  the  east  is  a 
region  of  inter-reef  sediments.  The  inner  shelf,  a  region 
of  recent  sedimentary  processes  (depositional) ,  contains  a 
variety  of  surficial  sediments  and  is  characterized  by  a 
gentle  slope  of  0.07°.  Farther  eastward  on  the  inner  shelf 
is  a  region  of  f oraminiferal  sands  with  sandwaves  common. 
Westward  of  Howell  Hook  the  slope  is  approximately  0.2°  and 
increases  through  a  series  of  steps  until  the  escarpment  is 
reached  with  slopes  as  much  as  45°  [Ref .  12]  . 

As  part  of  the  environmental  characterization  for 
Exercise  LWAD  99-1,  sediment  samples  were  collected. 
However,  the  results  were  not  available  in  time  for 
inclusion  in  this  study.  Preliminary  analysis  indicates 
that  Howell  Hook  is  overlain  by  a  layer  of  silty  sand,  while 
the  inter-reef  region  is  composed  of  a  clay  silt  mixture. 
East  of  the  inter-reef  region  the  sediment  is  composed  of 
fine  to  coarse  grain  sand  particles  [Ref.  13] . 

Based  on  the  geo-acoustic  modeling  methods  of  Hamilton 
and  Bachman  [Ref.  14],  a  characterization  of  the  reflective 
properties  of  the  sediments  can  be  accomplished.  From  Table 
1,  the  mean  grain  size  can  be  determined  from  a  description 


of  the  sediment  type.  Grain  size  can  then  be  correlated  to 
porosity  and  saturated  bulk  density  as  illustrated  in 
Figures  6  and  7 .  Porosity,  which  is  a  measure  of  the  amount 
of  voids  occupied  by  water  in  a  material,  can  then  be  used 
to  determine  the  compressional  wave  speed  of  the  sediment 
(Figure  8).  The  characteristic  impedance,  the  product  of 
the  sediment  density  and  wave  speed,  is  an  indicator  of  the 
reflective  nature  of  a  particular  sediment  composition. 

The  geo-acoustic  method  described  above  indicates  that 
the  degree  of  reflectivity  can  be  determined  from  the 
sediment  type.  In  particular,  reflectivity  is  proportional 
to  grain  size.  It  is  anticipated,  therefore,  that  the 
region  on  the  eastern  side  of  site  1  should  be  the  most 
reflective  and  should  have  higher  reverberation  levels.  The 
Howell  Hook  region  should  be  less  reflective,  but  still  more 
reflective  than  the  inter-reef  region  which  should  be  the 
least  reflective  and  have  correspondingly  lower 
reverberation  levels. 

2 .    Event  Weather  and  Oceanography 

All  acoustically  relevant  meteorological  and 
oceanographic  parameters  were  essentially  invariant  during 
the  12 -hour  duration  of  event  31,  making  for  near- laboratory 
conditions.  The  sky  was  free  of  clouds  with  wind  speeds 
less  than  5  m/s.  The  sea  surface  was  calm  with  a 
significant  wave  height  of  0.5-0.7  m.   The  surface  current 


measured  during  the  event  averaged  around  0.25  m/s  in  a 
northerly  direction.  [Ref.  15] 

A  total  of  23  CTD  observations  were  taken  throughout 
site  1  between  6-9  February.  Figure  9  shows  the  sound 
speed  profiles  (SSP)  calculated  from  the  CTD  observations. 
Little  spatial  or  temporal  variability  between  the  profiles 
is  observed  (<  1  m/s) ,  reflective  of  the  near  constant 
oceanic-atmospheric  forcing.  The  mixed  layer  is  seen  to 
oscillate  between  18  -  25  m.  The  downward  refracting  nature 
of  the  SSP  ensures  strong  bottom  interactions  of  the 
acoustic  signal.  [Ref.  15] 

3 .    Bottom  Backscattering  Measurements 

Bottom  backscattering  measurements  taken  with  a 
vertical  array  during  LWAD  99-1  [Ref.  16]  indicated  that  the 
bottom   scattering   strength    (Sb)    exhibited   a   sin (6) 

dependence,  not  the  sin2  (8)  dependence  of  Lambert's  Law. 
Previous  results  obtained  from  a  shallow  water  area  south  of 
Long  Island  [Ref.  5]  established  a  sin(0)  bottom  scattering 

strength  relationship  with  a  silty  clay  bottom  having  a 
surface  sediment  compressional  wave  speed  of  less  than  1600 
m/s.   Over  sandy  bottoms  Sb  suggested  a  sin2  (9)  dependence  in 

agreement  with  Lambert's  Law. 

It  is  common  practice  to  make  backscattering 
measurements  in  just  one  location  and  then  extend  those 
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results  to  a  wider  area  [Ref.  12],  but  this  practice  can  be 
misleading.  The  spatial  distribution  of  sediments  in 
shallow  water  is  generally  highly  variable  over  short 
spatial  scale  lengths.  Scanlon  et  al .  [Ref.  5]  showed  that 
the  spatial  scale  can  grade  from  hard  sand  to  mud/ silt  clay 
sediments  over  distances  less  than  500  m. 

In  this  thesis,  a  bottom  sediment  map  is  constructed 
based  upon  spatially-varying  changes  in  the  beam 
reverberation  data  from  an  AN/SQS-53C  tactical  active  sonar. 
While  the  results  are  considered  preliminary  due  to  a  lack 
of  high  resolution  bottom  samples,  they  do  indicate  a  high 
degree  of  spatial  variability  of  sediment  types  that,  only 
on  the  average,  agree  with  the  pre-exercise  wide-area 
sediment  assessment  [Ref.  12].  It  will  be  instructive  to 
compare  the  results  of  the  frequent  sediment  grab  samples 
taken  during  LWAD  99-1  with  the  results  of  this  thesis.  A 
conclusion  that  can  be  stated  now  is: 

Single  point  bottom  backscattering  measurements  can  not 
predict  the  high  spatial  sediment  variability  and  inverse 
techniques  must  be  developed  to  estimate  this  variability. 

Another  limitation  of  bottom  backscattering  analysis  is 
the  assumption  that  the  scattering  occurs  at  the  sediment- 
water  interface.  Recent  studies  [Ref.  17]  show  that  the 
scattering  occurs  from  inhomogeneities  in  the  sediment  sub- 
bottom  for  slow  (i.e.,  clay,  mud)  sediments  at  low  grazing 
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angles  (Figure  10)  .  Therefore,  a  full  wave  model  with  an 
accurate  three-dimensional  model  of  the  sub-bottom 
geoacoustic  properties  must  be  used  to  estimate  TL  in 
shallow  water.  Just  a  few  years  ago  it  did  not  seem 
possible  to  run  a  full  wave  model  at  3  500  Hz,  but  with 
today's  computing  power,  it  is  feasible.  Within  a  few  years 
it  will  be  commonplace  to  run  these  full  wave  models  in  near 
real  time  and  the  real  time  IT  sediment -mapping  algorithm 
developed  here  is  intended  for  this  future  application. 
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Bulk 

grain 

Sediment 

No. 

Mean 

grain  size 

Sand, 

Silt, 

Clay, 

density 

type 

samples 

mm 

phi  size 

% 

% 

% 

g/cm3 

Sand 

Coarse 

2 

0.5285 

0.92 

100.0 

0.0 

0.0 

2.710 

Fine 

28 

0.1638 

2.61 

92.2 

4.1 

3.7 

2.709 

Very  Fine 

16 

0.0988 

3.34 

81.0 

12.5 

6.5 

2.680 

Silty  Sand 

40 

0.0529 

4.24 

57.0 

30.9 

12.1 

2.677 

Sandy  Silt 

47 

0.0340 

4.88 

30.3 

57.8 

11.9 

2.664 

Silt 

19 

0.0237 

5.40 

7.8 

80.1 

12.1 

2.661 

Sand-silt-clay 

29 

0.0177 

5.82 

31.7 

42.9 

25.4 

2.689 

Clayey  silt 

105 

0.0071 

7.13 

7.4 

58.3 

34.3 

2.656 

Silty  clay 

54 

0.0022 

8.80 

3.9 

34.8 

61.3 

2.715 

Table  1.  Continental  terrace  (shelf  and  slope)  environment; 
Average  sediment  size  analyses  and  bulk  grain  densities. 
From  Ref .  [14] . 
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II.   TREATMENT  OF  DATA 

A.    DESCRIPTION  OF  LWAD  99-1  REVERBERATION  DATA 

The  beam  reverberation  data  used  in  this  study  was 
collected  using  the  AN/SQS-53C  fleet  tactical  sonar  aboard 
the  USS  Vicksburg,  during  event  31  of  exercise  LWAD  99-1  on 
9  February  1999  [Ref.  10].  Beam  reverberation  data  were 
analyzed  from  two  modes  of  operation  of  the  AN/SQS-53C,  the 
SD  (surface  duct)  and  VD  (variable  depression)  modes.  For 
consistency,  all  pings  that  were  analyzed  used  either  a 
linear  phase  modulated  pulse  or  a  CP  (coded  pulse)  waveform. 

The  SD-mode  uses  a  3  00  millisecond  pulse  length  and  has 
a  3  60°  azimuthal  transmit  coverage  employing  72  receive 
beams.  The  VD-mode  uses  a  500  millisecond  pulse  length  and 
has  a  12  0°  azimuthal  transmit  coverage  employing  24  receive 

beams.  In  both  modes,  the  3-dB  beamwidth  is  5  degrees.  The 
frequency  range  is  from  2400  to  4500  Hz.  The  depression 
angle  was  3°,  while  the  range  scale  was  set  to  20  kyds . 

Event  31  was  conducted  such  that  the  VD-mode  pings  were 
directed  toward  the  center  of  the  CG  box.  This  was  done  to 
get  complete  coverage  of  the  box  and  have  RL  information 
about  the  same  geographic  location  from  different 
directions.   When  the  CG  was  traversing  the  western  portion 

of  the  box,  on  a  course  of  000°T,  the  sonar  center  section 
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true  bearing  was  090°T  (Figure  11) ;  during  the  northern  leg 
(course  090°T)  the  sonar  center  section  true  bearing  was 
180°T  (Figure  12),  etc.   The  center  section  true  bearing  for 

the  VD-mode  is  the  bearing  between  beams  12  and  13 . 

All  SD-mode  pings  (omni-directional)  had  a  center 
section  true  bearing  of  180°T.  The  center  section  true 
bearing  for  the  SD-mode  is  the  bearing  between  beams  3  6  and 
37.   Beam  1  would  therefore  cover  from  000°T  to  005°T,  beam 

2  from  005°T  to  010°T,  etc.   Beam  72  would  cover  bearings 
355°T  to  360°T. 

B.    HYPOTHESIS 

Reverberation  is  a  combination  of  backscattered  energy 
from  within  the  water  column  (volume) ,  the  sea  surface  and 
the  sea  floor.  For  this  study  it  is  assumed  that  the  volume 
and  sea  surface  components  of  reverberation  were  constant. 
Therefore,  it  is  hypothesized  that  for  a  single  ping  from 
the  AN/SQS-53C  different  RL  beam  (bearing)  time  series 
should  only  differ  because  of  variations  in  backscattered 
energy  from  differing  bottom  types.  This  hypothesis  is 
further  extended  to  encompass  a  set  of  pings  (as  many  as  41 
pings  were  used  in  one  data  set)  ,  as  long  as  they  occur  in 
the  same  general  area  and  there  are  no  changes  to  the  ping 
type. 
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The  following  assumptions  were  made  to  support  this 
hypothesis : 

(1)  The  sound  speed  profile  in  the  area  around  the 
source  (within  -20  kyds)  is  constant. 

•  As  seen  in  Figure  9,  there  is  little  spatial  or 
temporal  variability  of  the  sound  speed  profile. 

(2)  The  output  source  level  (SL)  from  the  AN/SQS-53C  in 
any  direction  does  not  vary  (i.e.,  the  AN/SQS-53C  SL  is 
independent  in  azimuth) . 

•  The  difference  in  performance  of  the  AN/SQS-53C  in 
any  direction  is  considered  negligible  for  this 
experiment.   Analysis  of  individual  RL  time  series, 

though,  indicates  that  within  45°  of  the  CG  stern  an 

increase  in  ambient  (background)  noise  (AN)  due  to 
own  ship  noise  is  observed  (Figures  13  and  14)  . 
This  phenomenon  was  accounted  for  in  the  processing 
of  the  RL  data . 

(3)  There  are  no  major  topographic  changes  that  would 
significantly  change  the  RL  curve  (i.e.,  no 
upslope/downslope  effects,  no  pinnacles,  etc.) 

•  Based  on  the  geologic  characterization,  there  were 
no  bottom  topography  features  that  would 
significantly  effect  the  received  reverberation 
levels . 
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C.    PROCESSING  METHOD 

An  example  of  how  the  sediment  effects  reverberation 
can  be  seen  by  comparing  two  different  RL  time  series  in 
which  distinctly  different  characteristics  can  be  seen. 
Figure  13  is  a  single  beam  (beam  17)  RL  time  series  bearing 
082. 5°T.     One   sees   that   the  RL  curve  reduces   to  the 

background  noise  level  at  around  18  s  where  the  AN  level  at 
3.5  kHz  is  approximately  3  5  dB.  In  contrast,  Figure  15 
(beam  37,  bearing  182. 5°T)  shows  an  example  with  the  same  AN 

level,  for  which  the  RL  curve  does  not  reach  the  AN  level 
until  about  23.5  s.  Clearly  the  beam  data  is  RL-limited 
longer  in  Figure  15.  It  is  suggested  that  this  heightened 
RL  is  due  to  the  presence  of  a  more  reflective  bottom.  The 
relationship  between  RL  and  AN  as  a  function  of  detection 
range  is  illustrated  in  Urick  [Ref.  3,  Figure  2.2]. 

The  general  method  proposed  for  RL  curve  analysis  is  as 
follows:  (1)  determine  a  characteristic  RL  curve  for  an  area 
determined  from  the  average  of  many  pings  and  (2)  compare 
individual  beam  data  to  determine  its  difference  from  the 
area-wide  mean. 

A  variety  of  statistical  parameters  derived  from  the 
beam  RL  time  series  can  be  examined  to  determine  if  the 
characteristics  of  the  RL  curve  bear  any  correlation  to 
sediment  type.    The  statistical  parameters   investigated 
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include  the  mean  slope,  standard  deviation,  skew,  and  excess 

kurtosis . 

Before    continuing,    some    background    about    the 

statistical  method  used  in  this  thesis  is  required.    The 

statistics  calculated  are  not  based  on  the  entire  time 

series   but   rather   from   segments   of   the   time   series 

partitioned  into   0.5,   1,   or  2   second  increments.     In 

addition,  a  degree  of  overlap  of  the  data  is  often  useful, 

e.g.,  0.25,  0.5,  0.75  (or  25/50/75%  of  the  interval).  The 

statistics  derived  from  the  overlapped  segments  are  referred 

to  as  "moving"  statistics   (moving  average,   moving  skew, 

etc.)  .   For  example,  a  25-second  time  series,  with  320  data 

points  per  second,   contains   8000  data  points.   With  an 

interval /overlap  combination  of  1/0.75  the  first  data  point 

of  the  moving  mean  would  be  the  mean  of  points  1  to  32  0,  the 

second  point  would  be  the  mean  of  points  81  to  400,  the 

third  would  be  the  mean  of  points  161  to  480,  etc.   Figure 

16  shows  an  individual  RL  time  series  (black)  overlaid  by 

the  1/0.75  moving  mean  (white)  RL  curve. 

1.    Example  Using  Slope  Data,  SD-Mode,  1/0.75  Moving 
Statistics 

From  a  data  set  of  41  pings  an  area-wide  average  slope 

is  calculated  by  first  taking  the  1/0.75  moving  mean  for  all 

of  the  beams  (72)  in  each  transmitted  ping.   The  slope  of 

each  beam  is  determined  by  calculating  the  gradient  (d/dt) 

of  the  mean  RL  curve.   Calculating  the  mean  from  all  (72x41) 
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slope  data  sets  generates  the  area-wide  average  slope  data. 
The  area-wide  average  RL  slope  for  all  the  pings  examined  is 
comprised  of  over  2500  individual  RL  slope  time  series. 
Figure  17  shows  the  area-wide  average  (solid  line)  RL  slope 
and  its  standard  deviation  (--) . 

Figure  17  also  shows  the  beam  slope  for  a  particular 
beam/ping  relative  to  the  area-wide  average  slope.  Each 
point  along  the  single  beam  slope  has  an  associated  range 
and  bearing  which  corresponds  to  a  geographic  position.  By- 
determining  the  difference  between  the  area-wide  average 
slope  and  the  individual  beam  slope,  a  geographic  contour 
plot  can  be  generated  to  display  the  anomaly  from  the  mean 
(Figure  18)  .  Subsequently,  different  pings  can  be  compared 
to  determine  if  any  geographic  correlation  exists  between 
sets  of  pings  to  identify  a  unique  but  consistent  area 
exhibiting  similar  RL  characteristics. 

The  other  statistical  parameters  (STD,  skew,  and  excess 
kurtosis)  can  be  analyzed  using  an  identical  method  to  that 
used  for  the  slope. 
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III.  ANALYSIS 

A.    ANALYSIS  OF  RL  SLOPE  DATA 
1 .    SD-Mode 

A  thorough  analysis  of  all  the  RL  slope  plots  was 
conducted  to  create  a  geographic  representation  of  sediment 
type  based  on  consistent  spatial  variations  in  the  RL 
slopes.  From  a  continuous  sequence  of  41  SD-mode  pings  a 
contour  plot  was  created  based  on  the  deviations  of  the 
individual  RL  slopes  from  the  area-wide  average. 

Through  the  analysis  of  the  slope  data  and  comparisons 
to  that  from  different  pings,  RL  slope  features  can  be  seen 
that  are  consistent  in  their  geographic  location.  The  RL 
curves  from  ping  34  exhibit  slope  characteristics  that  are 
common  in  all  pings  in  the  set  and  will  be  discussed  in 
detail.  Figure  18  is  a  geographic  contour  plot  of  the 
deviations  of  each  of  the  72  azimuthal  beams  of  ping  3  4  RL 
slope  from  the  area-wide  average  RL  slope.  It  is 
anticipated  that  regions  where  the  deviations  (from  the 
average  RL  slope)  are  positive  correspond  to  more  reflective 
areas,  while  regions  with  a  negative  deviation  are  less 
reflective.  For  example,  Figure  19  shows  the  RL  and  slope 
data  of  ping  34,  beam  22.  At  approximately  8.5  s  a  region 
of  increased  RL  is  observed  with  a  corresponding  increase  in 
RL  slope.   Figure  20  depicts  the  region  of  beam  22  enlarged 
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to  show  the  RL  slope  details.  The  same  geographic  region  of 
increased  RL  can  be  seen  on  contour  plots  of  different  pings 
(Figures  21  -  23).  Note  that  the  RL  slope  anomaly  seen  here 
is  consistent  in  its  geographic  location  indicating  that  it 
is  due  to  a  geographic  feature,  and  not  to  a  random 
occurrence . 

Figure  24  is  similar  to  that  of  Figure  18  but  is 
derived  from  a  composite  plot  of  all  41  SD-mode  pings.  It 
is  the  average  deviation  of  RL  slope  from  the  area-wide 
average  RL  slope  at  each  geographic  location.  Being  an 
average  of  the  RL  slope  deviations,  this  plot  will  show  RL 
slope  features  that  have  geographic  consistency,  while 
eliminating   those   features   that   are  not   geographically 

consistent.   Note  that  the  area  within  +/-  45°  of  the  CG 

stern  were  not  included  in  the  composite  map  and  accounts 
for  the  notches,  especially  to  the  west,  seen  in  the  figure. 
An  analysis  of  Figure  24  shows  clear  similarities  to 
the  pre-exercise  sedimentary  characterization  (Figure  5)  . 
There  are  3  distinct  geographic  locations  observed  in  Figure 
24;  east  of  84.05°W  is  the  f oraminif eral  sands  area,  Howell 

Hook  is  west  of  84.3°W  and  the  middle  is  the  inter-reef  area 

which  is  interspersed  with  low  and  high  reflective  zones. 
The  highly  reflective  zones  could  be  uncovered  (or  thinly 
covered)  reef  segments  whereas  the  less  reflective  spots  may 
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be  due  to  depressions/basins  filled  with  unconsolidated 
sediments.  It  is  notable  that  the  foramini feral  sands  and 
Howell  Hook  regions  are  relatively  uniform,  while  the  inter- 
reef  region  shows  high  spatial  variability  of  sediment 
types . 

2 .    VD-Mode 

Due  to  the  limited  data  available  the  VD-mode  data  were 
not  used  to  create  a  composite  map  of  the  exercise  area. 
However,  individual  pings  were  analyzed  to  see  if  there  were 
any  geographically  consistent  RL  slope  features  as 
determined  from  the  analysis  of  the  SD-mode  data. 

Due  to  the  nature  of  the  exercise  the  VD-mode  pings 
were  directed  towards  the  center  of  the  CG  box.  This  was 
done  to  obtain  full  coverage  of  the  area  within  the  box  but 
unfortunately  created  limitations  for  this  study.  In 
particular,  sedimentary  features  such  as  Howell  Hook  and  the 
foramini feral  sands  region,  which  were  beyond  the  boundaries 
of  the  exercise  box,  would  not  visible  on  the  RL  plots. 

The  VD-mode  pings  did,  however,  demonstrate  similar 
geographically  consistent  features  as  seen  in  the  SD-mode. 
Figures  2  5  and  2  6  are  individual  VD-mode  pings  from  which 
geographic  features  can  be  seen.  Note  that  the  high  and  low 
bottom  reflection  areas  observed  in  Figures  2  5  and  2  6  are 
also  seen  at  the  same  location  in  the  SD-mode  composite  plot 
(Figure  24) . 
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B.    HOS  ANALYSIS 

The  same  analysis  methods  used  on  the  RL  slope  data 
were  employed  to  identify  the  correlation  of  HOS  to  sediment 
type.  Figures  27-29  are  plots  of  RL  and  the  corresponding 
STD,  skew,  and  excess  kurtosis.  The  RL  STD  and  skew  plots 
show  deviations  that  can  be  correlated  to  RL  features .  The 
RL  kurtosis  deviations,  however,  appeared  random  and  showed 
no  apparent  correlation  to  sediment  type.  The  observed  RL 
STD  and  skew  features  are  marginally  useful  in  that  they 
only  provide  an  inference  that  some  sedimentary  change  is 
present.  They  do  not  provide  information  about  the 
characteristics  of  that  feature  (i.e.,  whether  or  not  it  is 
due  to  a  more  or  less  reflective  sediment) . 

In  order  to  further  investigate  RL  statistics,  a 
distribution  analysis  was  conducted.  Figure  3  0  is  a  plot  of 
skew  versus  excess  kurtosis  for  the  detrended  RL  beam  data 
of  ping  189.  The  skew  and  excess  kurtosis  data  shown  are 
highly  correlated  with  negative  skew  (characterized  by  an 
extended  low  amplitude  tail  compared  to  a  Gaussian  curve) 
correlated  to  positive  excess  kurtosis  (characterized  by  a 
peaked  distribution  with  more  data  centered  around  the  mean 
compared  to  a  Gaussian  distribution) .  Preliminary 
distribution  analysis  indicates  that  the  reverberation  data 
is  not  only  non-Gaussian,  but  the  data  spans  several 
distribution  types  [Ref.  18]. 
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It  is  important  to  note  that  statistical  properties  of 
the  detrended  beam  time  series  are  not  commonly  analyzed  in 
classical  detection  theory  applications.  Instead,  the 
output  of  a  replica  correlator  is  normally  the  input  to  a 
statistical  analysis,  and  normalizations,  such  as  whitened 
matched  filters,  etc.,  attempt  to  remove  the  amplitude 
decay,  or  slope,  of  the  beam  reverberation  time  series  from 
the  analysis  [Ref .  18]  .  In  this  thesis  we  have  taken  a 
quite  different  approach  by  using  the  beam  time  series  slope 
to  invert  bottom  properties  with  the  HOS  of  the  detrended 
time  series  as  derived  product.  If  the  thesis  objective 
were  focused  on  detection  performance,  we  could  have  done  a 
replica  correlation  analysis  of  the  residual,  or  detrended 
beam  time  series.  Instead,  we  are  simply  trying  to 
correlate  the  HOS  of  the  detrended  time  series  to  the 
environment.  Since  little  is  known  about  the  HOS  of  the 
detrended  beam  time  series,  we  are  just  reporting  that  they 
were  highly  non-Guassian  and  further  investigation 
concerning  their  distribution  function  is  needed.  Although 
the  HOS  did  not  correlate  well  with  sediment  type,  this  does 
not  imply  that  they  will  not  correlate  with  some  other 
environmental  parameter,  such  as  water  depth  or  sound  speed 
profile. 
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C.    DISCUSSION 

1.  Tactical  Implications  of  Sediment  Map 

The  high  variability  of  bottom  sediment  type  in  shallow 
water  exerts  a  strong  tactical  significance  because  it  can 
unexpectedly  affect,  often  times  adversely,  the  results  of 
an  acoustic  search  which  is  executed  under  the  assumption 
that  the  bottom  sediment  is  homogenous.  Figure  31(A)  is 
indicative  of  a  search  plan  that  was  executed  assuming  that 
detection  ranges  were  constant  (regardless  of  sediment 
type)  .  This  "mowing  the  lawn"  technique  could  result  in  a 
lower  probability  of  detection  (PD)  than  anticipated  due  to 
regions  where  detection  ranges  are  degraded  because  of  a 
less  reflective  bottom.  Prior  to  an  acoustic  search  an 
environmental  update  could  be  conducted  which  includes  an  RL 
survey  from  which  a  sediment  map,  similar  to  the  one  in 
Figure  24,  could  be  used  to  modify  the  search  plan. 
"Searching  by  sediments",  as  seen  in  Figure  31(B),  would 
account  for  detection  range  differences  due  to  the  bottom 
sediment  type.  In  this  case,  the  acoustic  search  would  be 
executed  under  the  appropriate,  spatially-varying  PD. 

2 .  Tactical  Significance  of  Non-Gaussian  HOS 

The  analysis  of  the  detrended  beam  reverberation  data 
indicated  that  the  backscattered  energy  is  non-Gaussian, 
i.e.,  the  values  of  skew  and  excess  kurtosis  are  non-zero, 
and  this  also  has  important  tactical  implications.    The 
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signal  processing  design  incorporated  into  the  AN/SQS-53C 
and  other  sonar  systems  is  based  on  the  assumption  that  the 
RL  distribution  is  Gaussian.  Challenging  the  Gaussian 
assumption  raises  two  important  questions: 

(1)  Is   the   non-Gaussian   nature   of   the   beam 
reverberation   HOS   due   to   environmental   effects   (e.g., 
shallow  water  acoustic  propagation)   or  because  of  sonar 
signal  processing  design? 

(2)  If  the  RL  HOS  are  due  to  the  environment,  will  the 
HOS  change  when  calculated  for  a  different  environment 
(e.g.,   deep   water,   continental   slope   or   soft,   muddy 

bottoms) ? 

The  first  question  could  be  answered  by  performing  the 
same  statistical  analysis  of  AN/SQS-53C  RL  data  but  taken 
from  another  environment.  If  the  HOS  are  similar  in  varied 
environments,  the  nature  of  the  HOS  is  probably  due  to  the 
signal  processing  equipment.  In  this  case,  the  non-Gaussian 
issue  would  be  resolved  by  making  equipment  modifications. 

If  the  non-zero  HOS  are  in  fact  caused  by  environmental 
factors,  then  the  existing  signal  processing  algorithms, 
from  which  sonar  system  design  is  based,  are  not  valid. 
Currently  correlation  detectors  or  . energy  threshold 
detectors  are  optimized,  in  the  probability  of 
detection/probability  of  false  alarm  (PD/PFA)  sense  for 
Gaussian  beam  reverberation  statistics.    If  the  signal  is 
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not  Gaussian,  then  non-Gaussian  treatment  of  the  RL  data 
must  be  used  to  process  the  reverberation.  Solinsky  [Refs. 
19-21]  has  shown  that  non-zero  HOS  algorithms  can  be 
designed  which  will  efficiently  mitigate  the  effects  of 
reverberation . 

A  fundamental  conclusion  from  this  study  is  that  the 
statistical  distribution  and  its  associated  properties  of 
the  reverberation  HOS  needs  to  be  determined  in-situ  and 
incorporated  into  the  PD/PFA  algorithms  to  improve  echo  to 
reverberation  processing  of  the  AN/SQS-53C  sonar  system. 
The  HOS  characteristics  of  the  environment  can  be  quantified 
through  analysis  of  the  abundant  AN/SQS-53C  data  available 
in  the  Naval  Undersea  Warfare  Center  (NUWC)  active  sonar 
data  base. 

The  non-Gaussian  reverberation  statistics  problem  can 
be  addressed  in  several  ways,  including  the  development  of  a 
neural  net  (NN)  algorithm,  schematically  illustrated  in 
Figure  32.  The  NN  approach  can  be  adapted  to  the  changing 
environment  and  this  technique  is  highly  recommended  for 
further  development  for  the  AN/SQS-53C. 

Another  deterministic  approach  to  addressing  non- 
Gaussian  statistics  is  to  apply  matched  field  processing 
(MFP)  concepts.  Here  the  transmitted  signal  replica  is 
propagated  through  the  shallow  water  environment  using  an 
accurate  propagation  model  based  upon  accurate  environmental 
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parameters  determined  by  ITs  such  as  the  one  developed  in 
this  thesis.  Any  procedure  to  process  the  transmitted 
waveform  to  accurately  reflect  the  non-Gaussian 
characteristics  of  the  environment  will  improve  the 
correlation  of  the  received  energy  which  has  been  propagated 
in  the  environment  and  reflected  off  the  target  and  other 
scatterers.  If  the  transmitted  waveform  is  not  similar  to 
that  of  the  received  energy,  the  correlation  will  be  low  and 
detection  performance  poor.  On  the  other  hand,  if  the 
volume,  bottom  and  surface  scattered  returns  correlate  "too 
well"  with  the  transmitted  waveform,  the  active  sonar  system 
will  be  overloaded  with  false  targets.  MFP  or  neural  net 
algorithms  are  ways  to  address  the  problem  of  non-Gaussian 
statistics  by  including  the  realities  of  propagation  and 
scattering  in  shallow  water  in  the  processing  of  the 
transmitted  waveform  prior  to  correlation  with  received  beam 
data. 
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IV.   CONCLUSIONS  AND   RECOMMENDATIONS 

Exercise  LWAD  99-1,  which  took  place  in  February  1999, 
included  an  inversion  technique  experiment  during  which 
AN/SQS-53C  reverberation  data  was  collected.  The 
reverberation  data  was  used  to  investigate  the  possibility 
of  employing  inversion  techniques  to  infer  bottom 
sedimentary  properties. 

The  hypothesis  presented  in  this  thesis  was  based  on 
the  assumption  that  the  magnitude  of  backscattered  energy 
and  its  rate  of  spatial  decay  were  directly  related  to  the 
acoustic  reflectivity  of  the  seabed.  By  examining  the  slope 
of  the  RL  curves  from  various  sedimentary  provinces,  one 
should  be  able  to  identify  and  associate  various  reverberant 
areas  with  unique  sediment  types.  This  study  examined 
various  reverberation  level  statistics  including  RL  slope, 
standard  deviation,  skew,  and  excess  kurtosis,  to  test  this 
hypothesis.  Of  the  various  statistical  parameters  the  RL 
slope  was  found  to  be  most  clearly  aligned  with  the  degree 
of  bottom  reflectivity.  Over  less  reflective  (mud/silt) 
bottoms  a  decrease  in  RL  slope  was  noted  while  acoustic 
interactions  with  a  more  reflective  bottom  (sand)  causes  an 
increase  in  RL  slope.  Based  on  this  finding,  a  technique 
was  developed  using  the  deviation  of  the  RL  slope  for  a 
particular  ping/beam  from  an  area-wide  average  RL  slope  to 
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generate  a  geographic  map  illustrating  bottom  sedimentary 
characteristics . 

Analysis  of  higher  order  statistics  (HOS) ,  which 
included  standard  deviation,  skew  and  excess  kurtosis, 
revealed  no  clear  correlation  to  sediment  type.  Variations 
in  standard  deviation  and  skew  did  provide  indications  of 
sedimentary  change . 

An  analysis  of  the  statistical  distribution  of  the  HOS 
indicated  that  the  RL  data  did  not  derive  from  a  Gaussian 
population.  It  also  showed  that  the  distribution  did  not 
appear  to  come  from  any  well-known  distribution  (e.g., 
Poisson,  log  normal,  Rayleigh,  etc.).  The  non-Gaussian 
nature  of  the  RL  data  has  serious  implications  for  sonar 
signal  processing  because  most  signal  processing  algorithms 
assume  a  Gaussian  distribution  to  the  noise  field. 

It  is  recommended  that  the  RL  slope  technique  for 
generating  a  sediment  map  be  used  by  fleet  assets  as  part  of 
efforts  to  develop  in-situ  and  archived  bottom  reflection 
databases  and  that  these  results  be  incorporated  into 
shallow  water  ASW  search  planning.  In  addition,  further 
investigation  into  the  nature  of  RL  statistics  and  its 
impact  on  sonar  system  performance  should  be  undertaken. 
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APPENDIX  A.  FIGURES 


RL 


Combination  of  surface  and 
volume  reverberation 

Target  Echo 


Time/Range 
Surface  scattering 


Figure  1.  Schematic  and  characteristic  RL  curve  for  a 
deep  water  environment.   Reverberation  is  due  primarily 
to  backscattering  from  volume  and  sea  surface  scatterers 
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Combination  of  surface,  volume 
and  bottom  reverberation 


Target  echo  hidden  in  reverb 


Time/Range 


Surface  scattering 


SxbeEatron 


Shallow  Water 


Figure  2 .  Schematic  and  characteristic  RL  curve  for  a 
shallow  water  environment.   Reverberation  is  generally- 
louder  being  a  function  of  volume,  sea  surface,  and 
bottom  backscatterers .   High  RL  levels  often  mask  the 
target  echo. 
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Figure  3.   LWAD  99-1  Exercise  Area.   Event  31  (Inversion 
Techniques)  was  conducted  at  site  1.  Depths  in  meters. 
From  Ref .  [9]  . 
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Figure  4.   Event  31  experiment  geometry.   Over  900  pings 
of  the  AN/SQS-53C  tactical  sonar  were  completed  between 
0730Z  and  1800Z  on  9  February  1999.  From  Ref .  [10]. 


36 


LWAD  99-1  Event  31 


25.9 
25.85 

25.8 
25.75 

25.7 
1  25.65 

"(5 

_l 

25.6 
25.55 

25.5 
25.45 


!l 


-<__'_. 


25.4 
84 


?o 


In  ter  reef  sediment    L1 
P 
Foraminiferalsand 
with  sand  waves 

Central  Reef 
(Howell  Hook) 


84.4 


84.3 


84.2 
Longitude 


84.1 


83.9 


Figure  5.   Exercise  LWAD  99-1  surface  sediment 
distribution.  After  Ref.  [11]. 
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Figure  6.  Mean  grain  size  versus  porosity.  Round  dots 
represent  continental  terrace  (shelf  and  slope)  samples 
From  Ref .  [14] . 
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Figure   7.      Mean  grain   size  versus    saturated  bulk  density 
Round  dots   represent   continental    terrace    (shelf   and 
slope)    samples.    From  Ref .     [14]. 


39 


1900 
1860 
1820 
1780 


J2.  17^0 

E 


£i 


700 


LU 

>  1660 

Q 

Z 

O  1620 

CO 


1580 
1540  - 
1500  - 


•J  ' 

P  m 

•A 

•    t    • 


•  •  • 


%    V  •   •  • 
•       V 

•  •• 

•  X  •  < 

•  •  ••_#  • 


>•   •  •  /  • 


I       I      I 


1       I 


30        40  50         60         70  80         90 

POROSITY,  % 


100 


Figure  8.   Porosity  versus  sound  velocity,  continental 
terrace  (shelf  and  slope).  From  Ref.  [14]. 


40 


SITE  1  CTD  DOWNCASTS 
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Figure  9.  23  total  CTD  downcasts  were  taken  at  site  1 
from  2/6/99  to  2/9/99.  Bottom  interaction  of  acoustic 
energy  is  ensured  due  to  downward  refracting  nature  of 
sound  speed  profile  (SSP) .  From  Ref .  [15] . 
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Figure  10.   Scattering  mechanisms  are  from  (1)  the  water- 
sediment  interface,  (2)  sub-bottom  volume,  and  (3)  the 
sediment-basement  interface.  From  Ref.  [17]. 
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Figure  11.  VD-mode  ping  with  center  section  true  bearing 
of  090°T.  VD-mode  pings  are  steered  toward  the  center  of 
the  CG  box. 
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Figure  12 .  VD-mode  ping  with  center  section  true  bearing 
of  180°T.   VD-mode  pings  are  steered  toward  the  center  of 
the  CG  box. 
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Ping:  189,  Mode:  SD,  Beam:  17,  Brng(T):  82.5 
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Figure  13.   RL  time  series  for  a  forward  projected  beam. 
RL  reduces  to  the  ambient  noise  (AN)  level  of 
approximately  35  dB  at  about  18s.   Compare  with  Figure  14 
(stern  beam)  to  see  difference  in  AN  level. 
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Figure  14.   RL  time  series  of  stern  beam.   The  ambient 
noise  (AN)  level  is  approximately  2  0  dB  greater  than  non- 
stern  beam  due  to  ship  noise. 
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Ping:  189,  Mode:  SD,  Beam:  37,  Brng(T):  182.5 
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Figure  15 .  RL  curve  reduces  to  background  noise  level  at 
around  23.5  s.   Compared  to  figure  13  the  longer  RL  decay 
time  is  associated  with  backscattering  from  a  more 
reflective  bottom. 
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Ping:  189,  Mode:  SD,  Beam:  17,  Brng(T):  82.5,  Interval:  1 ,  Overlap:  0.75 
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Figure  16.   Individual  RL  time  series  (black)  overlaid  by 
the  1/0.75  moving  mean  (white)  RL  curve. 
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Ping:  189,  Mode:  SD,  Beam:  17,  Brng(T):  82.5,  Interval:  1 ,  Overlap:  0.75 
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Figure  17.   Slope  of  the  1/0.75  moving  mean  (-*)  for  an 
individual  beam  compared  with  the  area-wide  average  slope 
(solid  line)  and  its  standard  deviation  (--). 
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Figure  18.   Ping  34.   Geographic  contour  plot  of 
deviations  of  all  72  beams  of  slope  data  from  the  area- 
wide  average  slope.   Warmer  colors  represent  more 
reflective  sedimentary  areas  than  the  cooler  colors. 
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Ping:  034,  Mode:  SD,  Beam:  22,  Brng(T):  107.5,  Interval:  1,  Overlap:  0.75 
150, 1 i 1 
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Figure  19.   Ping  34,  beam  22.   An  increase  in  RL  at 
corresponds  with  an  increase  in  RL  slope. 
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Ping:0034,Mode:  SD,  Interval:  1,  Overlap:  0.75 
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Figure  20.   Ping  34,  beam  22.   Enlarged  contour  plot  to 
display  the  region  of  increased  RL. 
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Ping:  25  ,  Mode:  SD,  Interval:  1 ,  Overlap:  0.75 
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Figure  21.   Ping  25.   Note  region  of  increased  RL  at  same 
geographic  location  as  Ping  34. 
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Ping:0030,Mode:  SD,  Interval:  1 ,  Overlap:  0.75 
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Figure  22.   Ping  30.   Note  region  of  increased  RL  at  same 
geographic  location  as  Ping  34. 
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Ping:  35  ,  Mode:  SD,  Interval:  1,  Overlap: 0.75 
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Figure  23.   Ping  35.   Note  region  of  increased  RL  at  same 
geographic  location  as  Ping  34. 
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Figure  24.  Composite  plot  of  deviation  of  RL  slope  from 
the  area-wide  average  for  41  SD  mode  pings.   Warmer 
colors  represent  more  reflective  sedimentary  areas  than 
the  cooler  colors. 
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Figure  25.   VD-mode  ping  2.   Note  that  the  features 
pointed  out  are  also  seen  in  the  SD-mode  composite  plot 
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Figure  26.   VD-mode  ping  4.   Note  that  the  features 
pointed  out  are  also  seen  in  the  SD-mode  composite  plot 
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Ping:  034,  Mode:  SD,  Beam:  23,  BrngfT):  292.5,  Interval:  1 ,  Overlap:  0.75 
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Figure  27.  Ping  34,  beam  22.   An  increase  in  RL  at  8 . 5  s 
corresponds  to  an  increase  in  RL  STD. 
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Ping:  034,  Mode:  SD,  Beam:  23,  Brng(T):  292.5,  Interval:  1 ,  Overlap:  0.75 
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Figure  28.   Ping  34,  beam  22.   An  increase  in  RL  at  8.5  s 
corresponds  to  an  RL  skew  feature. 
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Ping:  034,  Mode:  SD,  Beam:  23,  Brng(T):  292.5,  Interval:  1 ,  Overlap:  0.75 
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Figure  29.   Ping  34,  beam  22.   The  RL  kurtosis  features 
appear  random  with  no  apparent  correlation  to  RL. 
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Figure  30.   Skew  vs.  excess  kurtosis  for  the  detrended 
beam  reverberation  data  of  ping  189. 
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Figure  31.   Old  search  method  (A)  vs.  new  search  method 
(B)  with  knowledge  of  local  sediment  types. 
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Figure  32.   Two-Dimension  feature  space.  From  Ref .  [21] 
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APPENDIX  B.  MATLAB  FILES 

The  matlab  program  files  enclosed  are  for  analyzing 
beamform  data  from  exercise  LWAD  99-1.  The  program  files  do 
not  include  the  actual  beamform  ping  data,  or  the  ping 
parameter  files.  The  matlab  files  included  here  are  for 
guidance  in  creating  a  graphical  user  interface  (GUI)  for 
analyzing  beamform  data. 

Beamform  data  files  contain  the  following  variables: 
bandwidth  -  bandwidth (Hz) 
bmf  -  complex  beamform  data 
centerfreq  -  center  frequency  (Hz) 
ping_date  -  date  (e.g.  '02/09.99') 
ping_time  -  time  (e.g.  '11:18:15') 
pulselength  -  pulse  length  (ms) 
samplerate  -  samplerate  (Hz) 

The  parameter  file  'evt31.mat'  contains  a  vector  of 
times  and  corresponding  CG  and  TGT  positions  from  LWAD  99-1 
event  31: 

evt31_time  -  time  in  seconds  from  midnight 
evt31_date  -  event  31  date  (e.g.  '02/09.99') 
cg_lat  -  vector  of  CG  latitudes 
cg_lon  -  vector  of  CG  longitudes 
tgt_lat  -  vector  of  TGT  latitudes 
tgt_lon  -  vector  of  TGT  longitudes 
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The  parameter  file  ' evt31prm.mat '  contains  a  vector  of 
times  and  corresponding  ping  parameters : 
botdepth  -  bottom  depth  (fm) 

course  -  CG  course  (°T) 

ctrbrng  -  Center  of  ping  section  bearing  (°T) 

deangle  -  sonar  depression  angle 

speed  -  CG  speed  (kts) 

timeday  -  time  in  seconds  from  midnight 
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%LWAD 

% 

%   Run  a  GUI  to  analyze  LWAD  99-1  data 

% 

%   A.  Begin  by  selecting  a  ping  to  analyze 

%   The  following  will  occur: 

%  (1)  A  map  of  the  LWAD  99-1  exercise  will  be  created 

%  showing  the  following: 

%  (a)  blue  eg  track  history,  red  circle  is  eg 

%  current  location 

%  (b)  red  tgt  track  history,  blue  circle  is  tgt 

%  current  location 

%  (c)  ping  section  is  shown  in  yellow  with  dashed 

%  yellow  center 

%  section  true  bearing  line.   Yellow  asterisks 

%  show  lOnm  range. 

%  (d)  latitude/ longitude  grid  lines  (black  dashed 

%  lines) 

%  (e)  Sediment  map  showing  Howell  hook  (black  circles) 

%  with  phi  size  approx  4-7,  clay  silt  (black  dots) 

%  with  phi  size  approx  6-8,  and  sand 

%  (black  slashed  lines)  with  phi  size  approx  1. 

%  This  map  is  based  on  data  received  from 

%  J.  Fulford  (NRL  Stennis) . 

% 

%  (2)  A  list  of  ping  parameters  will  be  displayed  in 

%  the  'LWAD  Parameters'  section. 

% 

%  (3)  A  data  analysis  section  will  be  created 

%  in  the  GUI. 

% 

%   B.  To  create  a  contour  plot  press  'Contour  Plot'  button. 

%  -  This  plot  is  created  by  determining  the 

%  difference  between  area  averaged  slope  data  and 

%  slope  of  individual  beam  data. 

%  -  The  plot  is  a  filled  contour  plot. 

%  -  A  regular (non-filled)  contour  plot  may  be 

%  created  by  using  the  function  ' lwad_contour .m' . 

% 

%    Note:  when  the  'contour  plot'  button  is  pressed,  the 

%  contour  plot  is  created 

%  using  an  interval  of  1  second  and  an  overlap  of 

%  .75.   The  interval /overlap  is  not  what  is  selected 

%  in  the  LWAD  Data  Processing  Program  GUI.   If  you 

%  want  to  investigate  different  interval /overlap 

%  combinations,  you'll  probably  want  to  run 

%  ' lwad_contourf '  function  separate  from  the  LWAD 
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%  Data  Processing  Program  GUI. 

% 

%  C .  To  analyze  beamf orm  data 

%  (1)  select  the  type  of  plot(s)  to  display 

%  -  201og|data|  -  is  RL  data  in  dB  with  averaged 

%  data  in  red 

%  -  Slope  -  slope  of  RL  data  in  blue  with  area 

%  averaged  slope  in  red.   Dashed  red 

%  line  is  +/-  1  standard  dev. 

%  -  STD  -  standard  deviation  in  blue  with  area 

%  averaged  std  in  red.   Dashed  red  line 

%  is  +/-  1  standard  dev. 

%  -  Skew  -  Skewness  in  blue  with  area  averaged 

%  skewness  in  red.   Dashed  red  line  is 

%  +/-  1  standard  dev. 

%  -  Kurtosis  -  Excess  Kurtosis  in  blue  with  area 

%  averaged  excess  kurtosis  in  red. 

%  Dashed  line  is  +/-  1  std. 

%  (2)  Select  the  beam  you  would  like  to  analyze 

%  1-24 (VD)  or  1-72 (SD) 

%  (3)  Select  the  interval  and  overlap  to  be  used  by 

%  mmean.m,  mstd.m,  mskewness.m,  mkurtosis.m.   I 

%  usually  use  1  sec  and  .75  overlap.   Other 

%  combinations  may  be  used  to  see  the  effect. 

%  (4)  Press  'Finish'  to  see  the  plot(s) . 

%  (5)  Press  'Close  plots'  to  close  all  plots  but  the 

%  event  map,  and  the  LWAD  GUI. 

% 

%  D.  To  analyze  another  file  press  'Browse'. 

% 

%  E.  Press  'Quit'  to  end  your  session  and  close  all 

%  open  plots . 

% 

%  See  also  LWAD_CONTOUR ,  LWAD_CONTOURF ,  LWAD_MAP, 

%  LWAD_POSIT,  LWAD_PRM,  MMEAN,  MSTD,  MKURTOSIS, 

%  MSKEWNESS . 

%  Created:  May  8,  1999 

%  David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%  David  Schalm 

menu_num  =  1 ; 

figure 

load  'e:\schalm\lwad  99-l\data_drive' 

lwad  menul 
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%LWAD_MENU1 

% 

%    This  will  create  a  UI  to  select  a  file  for  the  LWAD 

%    program 

% 

%    See  also  LWAD 

%  Created:  April  29,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

fig  =  gcf; 

set (fig, . . . 

'Units' , 'points' ,  ... 
'Color' ,[0.80.80.8],  ... 
' MenuBar ' , ' none ' ,  ... 
'NumberTitle' ,  'off,  ... 

'Name' , 'LWAD  Data  Processing  Program' ,  ... 
'Position' , [450  240  300  300],  ... 
'Resize'  ,  'off  )  ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment ' , 'left' ,  ... 
'Position' , [10  280  180  10],  ... 
'String' , 'LWAD  Data  file',  ... 
'FontWeight' , 'bold' ,  ... 
'Style' , 'text' ) ; 
hO  =  uicontrol ( 'Units ', 'points ' ,  ... 
'BackgroundColor' , [1  11],  ... 
'HorizontalAlignment ',' left ' ,  ... 
'Position' , [10  260  180  15],  ... 
' String' , [data_drive  '*.mat'],  ... 
'Style' , 'edit' ) ; 
callbackl  =  [' [data_f ile, data_drive]  =  ' . . . 

uigetfile(data_drive,  '  »"*.*'"  ');'  ... 
set (fig,'  '''Name'''  ','  ''' Loading  File ...'' '  .. 
) ; '  'if  data_file  ==  0,  load  '  ... 
' 'e:\schalm\lwad  99-l\data_drive' ' '  ';'... 
set (fig,'  '''Name'''  ','  ... 
' 'LWAD  Data  Processing  Program' ' '  ' ) ; '  ... 
,else,  lwad_menu2 ,  end']; 
uicontrol ( 'Units ', 'points ' ,  ... 
'Callback' , callbackl,  ... 
'Position' , [200  260  40  15] ,  ... 
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Ill 

I         I         III 


' String ' , ' Browse ' ) ; 
callbacks  =  [ 'ButtonName=questdlg ( '  ... 

Are  you  sure  you  want  to  quit?' ' ' . . . 

Quit  LWAD' ' '  ','  '''Yes'''  ','  ... 
' ' ' No ' ' '  ' , '  ' ' ' Yes ' ' '  ');'... 

'if  length (ButtonName)  ==  3,  close  all,  end']; 
uicontrol ( 'Units' , 'points' ,  ... 
'Callback' ,callback2,  ... 
'Position' , [10  10  40  15] ,  ... 
'String' , 'Quit' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.5  0.5  0.5],  ... 
'HorizontalAlignment ' , 'left' ,  ... 
'Position' , [5  250  290  1],  ... 
'String' , '  ' ,  ... 
'Style' , 'text' ) ; 
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%LWAD_MENU2 

% 

%  This  will  display  the  file  parameters,  and  add  a  UI 

%  section  for  selecting  plot(s)  and  plot  options. 

% 

%  Note:  when  the  'contour  plot'  button  is  pressed,  the 

%       contour  plot  is  created  using  an  interval  of  1 

%       second  and  an  overlap  of  .75.   The  interval/overlap 

%       is  not  what  is  selected  in  the  LWAD  Data  Processing 

%       Program  GUI.   If  you  want  to  investigate  different 

%       interval /overlap  combinations,  see 

%        ' lwad_contourf .m' .   You'll  probably  want  to  run 

%        ' lwad_contourf '  separate  from  the  LWAD  GUI. 

% 

%  See  also  LWAD,  LWAD_CONTOURF 

%  Created:  April  29,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

fig  =  gcf; 

save  'e:\schalm\lwad  99-l\data_drive .mat '  data_drive 
fname  =  [data_drive  data_f ile] ; 
load ( fname , ' -mat ' ) 

[course, speed, ctrTbrng,botdepth, deangle]  =  lwad_j?rm ( fname ) ; 
[cglat , cglon, tgtlat, tgtlon]  =  lwad_jposit (fname) ; 
num_of_beams  =  size (bmf ,  2 )  ; 
set (hO, ' string' , data_f ile) ; 

uicontrol ( 'Units ', 'points' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

'HorizontalAlignment '  ,  ' lef t ' ,  ... 

'Position' , [10  230  180  10],  ... 

'String' , 'LWAD  Parameters' ,  ... 

'FontWeight' ,  'bold',  ... 

'  Style '  ,  '  text '  )  ,- 
uicontrol ( 'Units' , 'points' ,  ... 

'BackgroundColor' , [0 .8  0.8  0.8],  ... 

'HorizontalAlignment ',' lef t' ,  ... 

'Position' , [10  210  140  10],  ... 

' String' ,[' Ping  date:  '  ping_date] ,  ... 

' Style ' , ' text ' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
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'HorizontalAlignment ' , ' lef t ' ,  ... 
'Position' , [150  210  140  10],  ... 
'String' ,[ 'Ping  Time:  '  ping_time; 
'Style' , 'text' ) ; 


lwad_map ( f name ) ; 

figure (fig) 

uicontrol ( ' Units ' , ' points 

' BackgroundColor ' , [0.8 

'HorizontalAlignment' , 

'Position' , [10  200  140 

'String' ,[ 'CG  Position 

num2str (cglon)  'W'],  ... 

'Style' , 'text' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8] 

'HorizontalAlignment ',' lef t' ,  . 

'Position' , [150  200  140  10],  .. 


0.8  0.8],  ... 
left',  ... 
10],  ... 
'  num2str (cglat) 


'N 


num2str (tgtlat)  'N 


' String' ,[ ' TGT  Position: 

num2str (tgtlon)  'W'],  ... 

'Style' , 'text' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' ,  [0.8  0.8  0.8],  ... 

' HorizontalAlignment ' ,  ' left '  ,  ... 

'Position' , [10  190  140  10],  ... 

'String' ,[ 'CG  Course (T) /Speed(kts) :  '  . 
num2str (course)  '/'  num2str (speed) ] , 

'Style' , 'text' ) ; 
uicontrol ( 'Units ',  'points '  ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

'HorizontalAlignment' , 'left' ,  ... 

'Position' , [150  190  140  10],  ... 

'String', ['Center  Bearing:  ' 

'Style' , 'text' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0 

'HorizontalAlignment' , 'left' 

'Position' , [10  180  140  10], 

'String', ['Bottom  Depth (fm) : 

'Style' , 'text' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0 

'HorizontalAlignment' , 'left' 

'Position' , [150  180  140  10], 

'String', ['Depression  Angle: 

'Style' , 'text' ) ; 
uicontrol ( 'Units' , 'points' ,  ... 


num2str (ctrTbrng) ] 


num2str (botdepth) ] , 


8] 


num2str (deangle) ] , 
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' BackgroundColor' , [0.8  0.8  0.8],  ... 

' HorizontalAlignment ' , ' left ' ,  ... 

'Position' , [10  170  140  10],  ... 

'String' ,[ 'Bandwidth:  '  num2str (bandwidth)  '  Hz'],  . 

'Style' , 'text' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0. 8  0.8  0.8],  ... 

'HorizontalAlignment ',' left ' ,  ... 

'Position' , [150  170  140  10],  ... 

' String' ,[ 'Center  Freq:  '  num2str (centerfreq)  ... 
'  Hz' ] , 'Style' , 'text' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

'HorizontalAlignment ',' left ' ,  ... 

'Position' , [10  160  140  10],  ... 

' String' ,[' Pulse  length:  '  num2str (pulselength)  ... 
'  msec ' ] , ' Style ' , ' text ' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

' HorizontalAlignment ' , ' left ' ,  ... 

'Position' , [150  160  140  10],  ... 

' String' ,[' Sample  rate:  '  num2str (samplerate)  ... 
'  Hz' ] , 'Style' , 'text' ) ; 

%  Create  time  vector 

data_start  =  1; 

data_duration  =  length(bmf) /samplerate  -  (1 /samplerate) 

+  data_start; 
time  =  (data_start : 1 / samplerate :data_duration) ' ; 

uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

' HorizontalAlignment ' , ' left ' ,  ... 

'Position' , [10  150  140  10],  ... 

' String' ,[ 'Data  Duration:  '  num2str (data_duration) . . 

'  Sec' ] ,  ... 
'Style' , 'text' ) ; 

uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.5  0.5  0.5],  ... 

' HorizontalAlignment ' , ' left ' ,  ... 

'Position' , [5  140  290  1],  ... 

'String' , '  ' ,  ... 

'Style' , 'text' ) ; 
uicontrol ( ' Units ' , ' points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 

' HorizontalAlignment ' , ' left ' ,  ... 
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'Position' , [10  120  180  10],  ... 
'String' 

'Select  the  plots  you  would  like  to  display. 
'FontWeight'  , 'bold',  ... 
'Style' , 'text' ) ; 
hi  =  uicontrol ( 'Units ', 'points' ,  ... 

'BackgroundColor' , [0 . 8  0.8  0.8],  ... 
'HorizontalAlignment ' , ' lef t ' ,  ... 
'Position' , [10  100  60  10] ,  ... 
'String', '20  log  |data|',  ... 
' Style ' ,  ' checkbox ' ) ; 
h2  =  uicontrol ( 'Units ', 'points' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment ',' lef t' ,  ... 
'Position' , [10  85  60  10] ,  ... 
'String' , 'Slope' ,  ... 
'Style',  'checkbox'); 
h3  =  uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment ',' lef t' ,  ... 
'Position' , [10  70  60  10],  ... 
'String' , 'STD' ,  ... 
' Style ' ,  ' checkbox ' ) ; 
h4  =  uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' ,   [0.8  0.8  0.8],  ... 
' HorizontalAlignment ' , ' left ' ,  ... 
'Position' ,  [10  40  60  10]  ,  ... 
' String ' , ' Kurtosis ' ,  ... 
'Style',  'checkbox'); 
h5  =  uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment' , 'left' ,  ... 
'Position' , [10  55  60  10],  ... 
'String' , 'Skew' ,  ... 
' Style ' ,  ' checkbox ' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
' HorizontalAlignment ' , ' left ' ,  ... 
'Position' , [85  100  65  10],  ... 
'String' , 'Beam  Selection: ' ,  ... 
'Style' , 'text' ) ; 
beam  =  {1:1: num_of _beams } ; 
h6  =  uicontrol ( 'Units ', 'points ' ,  ... 
'HorizontalAlignment ',' lef t' ,  ... 
'BackgroundColor' ,[111],  ... 
'Position' , [85  40  65  55],  ... 
' String' , beam,  ... 
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'Style' , 'listbox' ,  ... 
'ListBoxTop' ,  1) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment ' , ' lef t ' ,  ... 
'Position' , [160  100  65  10],  ... 
' String' ,' Interval  (sec)',  ... 
'Style',  'text'); 
interval_sel  =  {.5  1  2 } ; 
h7  =  uicontrol ( 'Units ', 'points ' ,  ... 
'BackgroundColor' ,[111],  ... 
'HorizontalAlignment ',' lef t' ,  ... 
'Position' , [160  40  60  55] ,  ... 
' String' , interval_sel,  ... 
'Style' ,  'listbox' ) ; 
uicontrol ( 'Units ', 'points ' ,  ... 

'BackgroundColor' , [0.8  0.8  0.8],  ... 
'HorizontalAlignment ',' lef t ' ,  ... 
'Position' , [230  100  65  10] ,  ... 
'String' , 'Overlap  (0-1)',  ... 
'Style' ,  'text' ) ; 
overlap_sel  =  {0  .25  .5  .75}; 
h8  =  uicontrol ( 'Units ', 'points ' ,  ... 
'BackgroundColor' ,[111],  ... 
'HorizontalAlignment ',' lef t' ,  ... 
'Position' , [230  40  60  55] ,  ... 
' String' , overlap_sel,  ... 
'Style',  'listbox'); 
callback4  =  ['set (fig,'  '''Name'''  ','  ... 
'  ' 'Getting  Contour  Plot.  . . '  "  ');'... 
' lwad_contourf ( f name , 1 , . 75 ) ; ' . . . 
'set (fig, '  ' ' 'Name' ' '  ' , '  ... 
' ' 'LWAD  Data  Processing  Program' ' '  ' ) ; ' ] 
uicontrol ( ' Units ' , ' points ' ,  ... 
'Callback' ,  callback4,  . . . 
'Position' , [80  10  60  15],  ... 
'String' ,  'Contour  Plot'); 
callback3  =  ['set (fig,'  '''Name'''  ','  ... 
' ' ' Getting  Plots ...'''  ');'... 
' lwad_menu3 ' ] ; 
uicontrol ( ' Units ' , ' points '  ,  ... 
'Callback' ,  callback3 ,  . . . 
'Position' , [250  10  40  15] ,  ... 
' String ' , ' Finish ' ) ; 
set (fig, 'Name' , 'LWAD  Data  Processing  Program') 
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%LWAD_MENU3 

% 

%   This  script  file  is  used  to  create  the  selected  plots 

% 

%    See  also  LWAD 

%  Created:  April  29,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

beam_s elect ion  =  beam{l} (get (h6, 'value' )) ; 
data_file  =  get (hO, ' string' ) ; 
interval  =  interval_sel { get (h7 , 'value' )} ; 
overlap  =  over lap_sel{ get (h8, 'value' )} ; 
spacing  =  interval  -  interval* overlap; 
N  =  round ( interval * samplerate ) ; 

%  Create  time  vector 

[mtime,ao]  =  mmean ( time, N, overlap) ; 

if  ao  ~=  overlap 

overlap  =  ao; 

disp(['Note:  actual  overlap  is  '  num2str (ao) ] ) ; 
end 

%%%%%%%   PLOT    THE    DATA    %%%%%%%%%%%%% 

%  Display  plots  of  data/mean,  slope,  kurtosis,  skewness 

if  data_file(5:6)  ==  'SD' 

area_s tat s_f name  =  ['e:\schalm\lwad  99-l\area_stats\ ' . . 
data_file(5 :6)  '_stats'  ... 
num2str (interval)  '.mat']; 
else 

area_s tat s_f name  =  ['e:\schalm\lwad  99-l\area_stats\ ' . . 
data_file(19:20)  '_stats'  ... 
num2str (interval)  '.mat']; 
end 

load ( area_s tats_f name ) 

if  get (h8, 'value' )  ==  1 

area_stats  =  stats_0; 
elseif  get (h8, 'value' )  ==  2 

area_stats  =  stats_2  5; 
elseif  get (h8, 'value' )  ==  3 

area_stats  =  stats_5; 
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elseif  get (h8, 'value' )  ==  4 

area_stats  =  stats_75; 
end 

log_data  =  20*logl0 (abs (bmf ( : ,beam_s elect ion) ) ) ; 
thetal  =  ctrTbrng  -  num_of_beams*5/2  +  5/2; 
theta  =  thetal  +  5* (beam_selection  -  1); 
if  theta  <  0,  theta  =  theta  +  3  60;  end 

title_string  =  ['Ping:  '  data_f ile (14 : 16)  ',  Mode:  SD'  ... 
',  Beam:  '  num2str (beam_selection) . . . 
',  Brng(T):  '  num2str (theta)  ... 

',  Interval:  '  num2str (interval)  ',  Overlap:  '  ... 
num2str (overlap) ] ; 

%  plot  201og|data| 

if  get (hi, 'value' )  ==  1 
figure,  b  =  gcf; 

plot (time, log_data, 'k' ) ,  title ( title_string) , 
ylabel ( 'Reverberation  Level  (dB) ' ) ,  xlabel ( 'Time  (sec)') 
axis([0  data_duration  0  150]),  grid  on,  hold  on 
[mean_data]  =  mmean(log_data,N, overlap) ; 
plot (mtime,mean_data, 'w' ) 

end 

%  plot  the  slope  of  201og|data| 
if  get (h2, 'value' )  ==  1 

figure,  b  =  gcf; 

[mean_data]  =  mmean(log_data,N, overlap) ; 

slope_data  =  gradient (mean_data) ; 

plot (mtime, slope_data/ spacing, 'k*- ' ) , title (title_string) 

ylabel ('RL  Slope '), xlabel (' Time  (sec)'), grid  on,  hold  on 

plot (mtime,  area_stats (1, : ) /spacing, 'k' ) 

plot (mtime,  (area_stats ( 1 , : ) -area_stats (2  ,:))/.. . 
spacing, ' --k' ) 

plot (mtime,  (area_stats ( 1 , : ) +area_stats (2 ,:))/.. . 
spacing, ' --k' ) 

legend ([' Beam  '  num2str (beam_selection) . . . 

'  Slope' ] , 'Average  Slope' , 'Average  Slope  STD' ) 
end 

%  de trend  the  data  w/  a  3rd  order  polynomial 

log_data  =  log_data-polyval (polyf it ( time, log_data, 3 ) , time) ; 

%  Plot  the  STD 
if  get (h3, 'value' )  ==  1 
figure,  b  =  gcf; 
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[std_data]  =  mstd(log_data,N, overlap) ; 
plot (mtime, std_data, 'k*-' ) ,  title (title_string) 
ylabel ( 'STD' ) ,  xlabel ( ' Time  (sec)') 
grid  on,  hold  on 

plot (mtime,  area_stats (3 , : ) , 'k' ) 

plot (mtime,  (area_stats (3 , : ) -area_stats (4 , : ) ) , ' --k' ) 
plot (mtime,  (area_stats (3 , : ) +area_stats (4 ,:)),' --k' ) 
legend ([ 'Beam  '  num2str (beam_selection)  '  STD'],... 
'Average  STD', 'STD  of  Average  STD') 
end 

%  Plot  the  Kurtosis 

if  get (h4, 'value' )  ==  1 
figure,  b  =  gcf; 

[kurtosis_data]  =  mkurtosis (log_data, 1,N, overlap) ; 
plot (mtime, kurtosis_data, ' k*- ' ) , title ( title_string) 
ylabel ( 'Kurtosis ') ,  xlabel ('Time  (sec)') 
grid  on,  hold  on 

plot (mtime,  area_stats (7 , : ) , 'k' ) 

plot (mtime,  (area_stats (7 , : ) -area_stats (8 , : ) ) , ' --k' ) 
plot (mtime,  (area_stats (7 , : ) +area_stats (8 , : ) ) , ' --k' ) 
legend ([' Beam  '  num2str (beam_selection)  '  Kurtosis'], 
'Average  Kurtosis ',' STD  of  Average  Kurtosis') 

end 

%  Plot  the  Skew 

if  get (h5, 'value' )  ==  1 

figure,  b  =  gcf; 

[skew_data]  =  mskewness (log_data,N, overlap) ; 

plot (mtime, skew_data, 'k*- ' ) , title (title_string) 

ylabel ( ' Skew ' ) , xlabel ( ' Time  ( sec ) ' ) 

grid  on,  hold  on 

plot (mtime,  area_stats (5 , : ) , 'k' ) 

plot (mtime,  (area_stats (5, : ) -area_stats (6 ,:)),' --k' ) 

plot (mtime,  (area_stats (5, : ) +area_stats (6 ,:)),' --k' ) 

legend ( [ 'Beam  '  num2str (beam_selection)  '  Skew' ] , . . . 
'Average  Skew', 'STD  of  Average  Skew') 
end 

set (fig, 'Name' , 'LWAD  Data  Processing  Program' ) 

figure (fig) 

uicontrol ( ' Units ' , ' points ' ,  ... 

'Callback',  [ 'closebut ( [f ig  100] )  ' ]  ,  ... 

'Position' , [165  10  60  15],  ... 

'String' , 'Close  Plots'); 
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function  [area_stats]  =  lwad_stats (file_type, interval .. . 

, overlap) 
• %LWAD_STATS   Determine  the  statistics  for  the  LWAD  99-1 
%    exercise  area. 

%    [AREA_STATS,ERROR_LOG]  =  LWAD_STATS (FILE_TYPE, 
%    INTERVAL , OVERLAP )  will  return  the  statistics  for  the 
%    selected  parameters.  The  function  will  calculate 
%    statistics  for  each  individual  beam,  average  these 
%    statistics,  and  determine  the  standard  deviation. 
% 

%    FILE_TYPE  is  either  'sd',  'cp',  or  'cw'  ...  depending 
%    on  which  mode  you  want  to  analyze . 
% 

%    INTERVAL  and  OVERLAP  are  required  by  funtions  MMEAN, 
%    MSKEWNESS,  and  MKURTOSIS  to  calculate  the  respective 
%    statistics.   If  no  INTERVAL  or  OVERLAP  are  specified, 
%    the  program  defaults  to  INTERVAL  =  1,  OVERLAP  =  .75. 
% 

%   The  mean  and  std  for  the  following  statistics  are 
%    determined: 
%        (1)  Slope 
%        (2)  STD 
%        (3)  Skewness 
%        (4)  Excess  Kurtosis 
% 

%   Note:  STD,  Skewness,  and  Excess  Kurtosis  are  determined 
%         from  'detrended'  data.   The  data  is  detrended  by 
%         subtracting  a  best  fit  3rd  order  polynomial  from 
%         the  data.   Slope  is  not  detrended. 
% 

The  output  file  AREA_STATS  is  organized  as  such: 

mean  of  the  slope 

std  of  the  slope 

mean  of  the  std 

std  of  the  std 

mean  of  the  skewness 

std  of  the  skewness 

mean  of  the  kurtosis 

std  of  the  kurtosis 


%        Row  1 

%        Row  2 

%        Row  3 

%        Row  4 

%        Row  5 

%        Row  6 

%        Row  7 

%        Row  8 

% 

%  The  number  of  columns  will  be  a  function  of  the 

%  selected  INTERVAL  and  OVERLAP. 

% 

%  If  the  time  series  data  contains  a  zero  value  this  will 

%  cause  errors  in  the  statistics  because  log(O)  =  -inf. 

%  All  beams  which  contain  a  zero  value  are  discarded. 

%  An  error  log  may  be  obtained  as  an  output. 
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%  Created:  May  25,  1999 

%         David  Schalm 

% 

%  Last  Modified:  May  25,  1999 

%         David  Schalm 

load  data_drive 

drive_num  =  data_drive (1) ; 

if  nargin  ==  0 

file_type  =  input (' Select  file  type  (sd, cp, cw) . . . ' ) ; 
end 

if  file_type  ==  'sd' 

f ilepath  =  { [drive_num  . . . 

' :\LWAD  99-1  Data\Tape  17\sd\' ]  ... 
[drive_num  ' : \LWAD  99-1  DataXTape  18\sd\']}; 
elseif  file_type  ==  'cp' 

f ilepath  =  { [drive_num  . . . 

':\LWAD  99-1  DataXTape  17\vd\cp\']  ... 
[drive_num  ' : \LWAD  99-1  DataXTape  18\vd\cp\ ' ] } ; 
elseif  file_type  ==  'cw' 

f ilepath  =  { [drive_num  . . . 

':\LWAD  99-1  DataXTape  17\vd\cw\']  ... 
[drive_num  ' : \LWAD  99-1  DataXTape  18\vd\cw\ ' ] } ; 
end 

row_ index  =  0 ; 
error_index  =  0 ; 

for  b  =  1:1:2 
clear  fname 

D  =  dir (char (f ilepath (b) )) ; 
[ fname {1 : length (D) } ]  =  deal (D. name) ; 

%  get  rid  of  ' . '  and  ' . . '  filenames 
fname  =  fname (3 : length (fname) ) ; 

for  m  =  1 : 1 : length (fname) 

load ( [char ( f ilepath (b) )  char ( fname (m) ) ] ) 

disp ([ 'Current  file  is  '  filepath{b}  fname{m}]) 

X  =  (1:1: length (bmf ) ) ' ; 

N  =  round (samplerate* interval) ; 

num_of_beams  =  size (bmf , 2) ; 

%  Create  time  vector 
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data_start  =  1; 

data_duration  =  length(bmf ) /samplerate  - 
(1/samplerate) . . . 

+  data_start; 
time  =  (data_start : 1/samplerate :data_duration) ' ; 
[mtime,ao]  =  mmean (time, N, overlap) ; 

for  n  =  1 : 1 :num_of_beams 

%  Skip  over  data  with  0  values  -->  cause  errors 

%  in  statistics 

if  sum(bmf(:,n)  ==  0)  >  0 

else 

row_index  =  row_index  +  1; 

%  Calculate  amplitude  of  data  (dB) 
log_data  =  2  0*logl0 (abs (bmf ( : , n) ) ) ; 

mean_data  =  mmean (log_data,N, overlap) ; 
lwad_s tat s_s lope (row_index, : )  =  ... 
gradient (mean_data) ; 

%  Detrend  the  data 
data  =  log_data  -  ... 

polyval (polyf it (X, log_data, 3 ) ,X) ; 

lwad_stats_std(row_index, : )  =  ... 

mstd (data, N, overlap) ; 
lwad_stats_skewness (row_index, : )  =  ... 

mskewness (data, N, overlap) ; 
lwad_stats_kurtosis (row_index, : )  =  ... 
mkurtosis (data, 1,N, overlap) ; 
end 
end 
end 
end 

[M,N]  =  size (lwad_s tat s_s lope) ; 

lwad_stats_slope ( (M+l) , : )  =  mean (lwad_stats_slope) ; 
lwad_s tat s_s lope ( (M+2) , : )  =  std(lwad_stats_slope (1 :M, :  )  )  ; 

lwad_stats_std( (M+l) , : )  =  mean (lwad_s tat s_std) ; 
lwad_stats_std( (M+2) , : )  =  std(lwad_stats_std(l :M, : ) ) ; 

lwad_stats_skewness ( (M+l) , : )  =  mean (lwad_stats_skewness) ; 
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lwad_stats_skewness ( (M+2) , : )  =  ... 
std(lwad_stats_skewness (1  :M,  : ) )  ; 

lwad_stats_kurtosis ( (M+l) , : )  =  mean (lwad_stats_kurtosis) ; 
lwad_stats_kurtosis ( (M+2) , : )  =  ... 
std(lwad_stats_kurtosis (1 :M, : ) ) ; 


area_stats (1:2, 
area_stats (3:4, 
area_stats (5:6, 
area_stats (7:8, 


)  =  lwad_stats_slope( (M+l) : (M+2) ,  :  )  ; 

)  =  lwad_stats_std( (M+l) : (M+2) , : ) ; 

)  =  lwad_stats_skewness ( (M+l) : (M+2) , : ) ; 

)  =  lwad_stats_kurtosis ( (M+l) : (M+2) ,  :  )  ; 
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function  [y, actual_overlap]  =  mmean (x, N, overlap) 

%MMEAN   Moving  mean  value . 

%  MMEAN (x,N, overlap)  should  be  used  to  calculate  a  "moving 

%  average"  where  N  is  the  number  of  points  to  use  for  each 

%  interval,  and  OVERLAP  is  the  ammount  of  overlap  to  use. 

%  In  some  cases  an  assigned  overlap  may  not  result  in  an 

%  integer  value  for  STEP  (i.e.  number  of  points  to  advance 

%  after  each  calculation) .   In  this  case,  OVERLAP  will  be 

%  rounded  to  the  nearest  value  which  will  result  in  an 

%  integer  value  for  STEP.   The  ACTUAL_OVERLAP  can  be 

%  displayed  if  desired. 

% 

%  See  also  MEAN,  MSTD,  MKURTOSIS,  MSKEWNESS,  LWAD. 

%  Created:  May  8,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  2 

overlap  =  0 ; 
end 

step  =  round ( (1-over lap) *N) ;  %  Number  of  points  to  advance 

i  =  1; 

while  (N+ (i-1) *step)  <=  length(x) 

y(i)  =  mean (x((l+ (i-1 ) *step) : (N+( i-1) *step) )) ; 

i  =  i  +  1  ; 
end 

if  nargout  ==  2 

actual_overlap  =  l-(step/N); 
end 
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function  [y, actual_overlap]  =  mstd(x,N, overlap) 
%MSTD  Moving  standard  deviation 
% 

%   MSTD(x,N, overlap)  should  be  used  to  calculate  a  "moving 
%    std"  where  N  is  the  number  of  points  to  use  for  each 
%    interval,  and  OVERLAP  is  the  ammount  of  overlap  to  use. 
%    In  some  cases  an  assigned  overlap  may  not  result  in  an 
%    integer  value  for  STEP  (i.e.  number  of  points  to  advance 
%    after  each  calculation) .   In  this  case,  OVERLAP  will  be 
%    rounded  to  the  nearest  value  which  will  result  in  an 
%    integer  value  for  STEP.   The  ACTUAL_OVERLAP  can  be 
%   displayed  if  desired. 


% 


See  also  MEAN,  MEDIAN,  STD,  MIN,  MAX,  COV,  LWAD. 


%  Created:  May  8,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  2 

overlap  =  0; 
end 

step  =  round ( (1-overlap) *N) ;  %  Number  of  points  to  advance 

i  =  1; 

while  (N+ (i-1) *step)  <=  length (x) 

y(i)  =   std(x( (l+(i-l)*step) : (N+ (i-1) *step) ) ) ; 

i  =  i  +  1  ; 
end 

if  nargout  ==  2 

actual_overlap  =  l-(step/N); 
end 
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function  [s, actual_overlap]  =  mskewness (x,N, overlap) 

%MSKEWNESS  Moving  Skewness . 

%   MSKEWNESS (x,N, overlap)  should  be  used  to  calculate  a 

%    "moving  skewness"  where  N  is  the  number  of  points  to  use 

%    for  each  interval,  and  OVERLAP  is  the  ammount  of  overlap 

%    to  use.   In  some  cases  an  assigned  overlap  may  not 

result 

%    in  an  integer  value  for  STEP  (i.e.  number  of  points  to 

%    advance  after  each  calculation) .   In  this  case,  OVERLAP 

%   will  be  rounded  to  the  nearest  value  which  will  result 

%    in  an  integer  value  for  STEP.   The  ACTUAL_OVERLAP  can  be 

%    displayed  if  desired. 


% 


See  also  SKEWNESS,  MMEAN,  MSTD,  MKURTOSIS,  LWAD, 


%  Created:  May  8,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  2 

overlap  =  0 ; 
end 

step  =  round ( (1-over lap) *N) ;  %  Number  of  points  to  advance 

i  =  1; 

while  (N+ (i-1) *step)  <=  length(x) 

s(i)  =  skewness (x( (1+ (i-1) *step) : (N+ (i-1) *step) )) ; 

i  =  i  +  1 ; 
end 

if  nargout  ==  2 

actual_overlap  =  l-(step/N); 
end 
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function  [k, actual_overlap]  =  mkurtosis (x, flag, N, overlap) 

%MKURTOSIS  Moving  Kurtosis . 

%  MKURTOSIS (x, flag, N, overlap)  should  be  used  to  calculate 

%  a  "moving  kurtosis"  where  N  is  the  number  of  points  to 

%  use  for  each  interval,  and  OVERLAP  is  the  ammount  of 

%  overlap  to  use.   In  some  cases  an  assigned  overlap  may 

%  not  result  in  an  integer  value  for  STEP  (i.e.  number  of 

%  points  to  advance  after  each  calculation) .   In  this 

%  case,  OVERLAP  will  be  rounded  to  the  nearest  value 

%  which  will  result  in  an  integer  value  for  STEP.   The 

%  ACTUAL_OVERLAP  can  be  displayed  if  desired. 

% 

%  MKURTOSIS (x, 1,N, overlap)  calculates  excess  kurtosis. 

% 

%  See  also  KURTOSIS,  MSTD,  MSKEWNESS,  MMEAN,  LWAD. 

%  Created:  May  8,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  3 

overlap  =  0 ; 
end 

step  =  round ( (1-overlap) *N) ;  %  Number  of  points  to  advance 

i  =  1; 

while  (N+ (i-1) *step)  <=  length(x) 
if  flag  ==  1 

k(i)  =  kurtosis(x( (l+(i-l) *step) : (N+ (i-1) *step) ) )  -  3 ; 
else 

k(i)  =  kurtosis (x( (1+ (i-1) *step) : (N+ (i-1) *step) )) ; 
end 

i  =  i  +  1; 
end 

if  nargout  ==  2 

actual_overlap  =  l-(step/N); 
end 
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function  [cglatout , cglonout , tgtlatout, tgtlonout]  =  ... 

lwad_jposit  (data_f  ile) 
%LWAD_POSIT 

%   Determine  the  coordinates  of  the  DD  (lat/lon)  based  on 
%    the  ping  time. 
% 
%    See  also  LWAD 

%  Created:  May  15,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

load (data_f ile, '-mat') 

load  evt31.mat 

ref_time  =  str2num(ping_time (1 : 2) ) *3600  +  ... 

str2num(ping_time (4 : 5) ) *60  +  str2num(ping_time (7 : 8) ) ; 

i  =  1; 

while  evt31_time (i)  <=  ref_time 
i  =  i  +  1; 

end 

interp_f actor  =  (evt31_time (i)  -  ref_time) /60; 
cglatout  =  cg_lat(i)  -  (cg_lat(i)  -  cg_lat (i-1) ) * . . . 

interp_f actor; 
cglonout  =  cg_lon(i)  -  (cg_lon(i)  -  cg_lon (i-1) ) * . . . 

interp_f actor ; 
tgtlatout  =  tgt_lat(i)  -  (tgt_lat(i)  -  tgt_lat (i-1) ) * . . . 

interp_f actor; 
tgtlonout  =  tgt_lon(i)  -  (tgt_lon(i)  -  tgt_lon(i-l) ) * . . . 

interp_factor; 
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function  [course, speed, ctrbrng, bo tdepth, deangle]  =  ... 

lwad_j?rm  (data_f  ile) 
%LWAD_PRM 

%    This  function  will  access  the  parameter  file 
%    ' evt31prm.mat '  to  determine  the  ping  parameters 
%    associated  with  the  selected  data  file. 
% 

%     course:   ships  course  in  degrees  True 
%      speed:   ships  speed  in  knots 

%    ctrbrng:   center  of  ping  section  in  degrees  True 
%    botdepth:   bottom  depth  in  fathoms 
%     deangle:   sonar  depression  angle  in  degrees 
% 
%    See  also  LWAD 

%  Created:  May  11,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

load (data_f ile, '-mat') 

load  evt31prm.mat 

ref_time  =  str2num(ping_time (1 :2) ) *3  600  +  ... 

str2num(ping_time (4 : 5) ) *60  +  str2num(ping_time (7 : 8) ) ; 

i  =  1; 

while  timeday(i)  <   ref_time 

i  =  i  +  1; 
end 
i  =  i-1; 

course  =  course (i); 
speed  =  speed (i) ; 
ctrbrng  =  ctrbrng (i); 
botdepth  =  botdepth ( i) ; 
deangle  =  deangle (i); 
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function  lwad_map ( f name ) 

%LWAD_MAP 

%  Create  a  map  of  the  exercise  area 

%  (1)  A  map  of  the  LWAD  99-1  exercise  will  be  created 

%  showing  the  following: 

%  (a)  blue  eg  track  history,  red  circle  is  eg 

%  current  location 

%  (b)  red  tgt  track  history,  blue  circle  is  tgt 

%  current  location 

%  (c)  ping  section  is  shown  in  yellow  with  dashed 

%  yellow  center  section  true  bearing  line. 

%  Yellow  asterisks  show  lOnm  range. 

%  (d)  latitude/ longitude  grid  lines  (black  dashed 

%  lines) 

%  (e)  Sediment  map  showing  hal  hook  (black  circles) 

%  with  phi  size  approx  4-7,  clay  silt  (black  dots) 

%  with  phi  size  approx  6-8,  and  sand 

%  (black  slashed  lines)  with  phi  size  approx  1. 

%  This  map  is  based  on  data  received  from 

%  J.  Fulford  (NRL  Stennis) . 

% 

%  If  no  filename  is  specified,  a  gui  will  open  to  allow 

%  the  user  to  select  a  file. 

% 

%  See  also  LWAD 

%  Created:  May  15,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  0 

load  data_drive 

[filename, pathname]  =  uigetf ile ( [data_drive] ) ; 

fname  =  [pathname  filename] ; 

data_drive  =  pathname; 

save  'e:\schalm\lwad  99-l\data_drive.mat '  data_drive 
end 

load ( fname , ' -mat ' ) 

[eglat , eglon, tgtlat, tgtlon]  =  lwad_posit (fname) ; 

[course, speed, ctrTbrng, bo tdepth, deangle]  =  lwad_j?rm ( fname ) ; 

load  evt31.mat; 

num_of_beams  =  size (bmf , 2 ) ; 


figure (100) 


89 


elf 

set (100, . . . 

'NumberTitle' ,  'off,  ... 

' Name ' , ' LWAD  99-1  Event  31'); 

%  Plot  the  sediment 

pic  =  imread( ' f lorida' , 'bmp' ) ; 

sed_lon  =  linspace (-84 . 5, -83 . 9 , 600) ; 

sed_lat  =  linspace(25.4,25.9, 600) ; 

image (sed_lon, sed_lat, pic) ; 

hold  on,  grid  on,  axis  xy 

%plot  eg  and  tgt  track  history- 
plot  (-cg_lon,  cg_lat)  ,  hold  on 

xlabel ( ' Longitude  (W) ' ) ,  ylabel ( ' Latitude  (N) ' ) 
title(['LWAD  99-1  Event  31']),  axis  equal,  axis  image 
plot (-eglon, eglat, 'ro' ) 
plot (-tgt_lon, tgt_lat, 'r' ) 
plot (-tgt Ion, tgtlat, 'bo' ) 

%  Plot  a  lOnm  range  circle,  and  ctr  bearing  line 
theta_one  =  ctrTbrng  -  num_of_beams*5/2 ; 
theta  =  theta_one: 5 : (theta_one+ (num_of_beams-l) *5) ; 
theta (num_of_beams+l)  =  theta_one; 

for  i  =  1 : 1 : (num_of_beams+l) 

circle_lat (i)  =  eglat  +  10*cos (theta (i) *pi/180) /60; 

circle_lon(i)  =  -(eglon  -  10*sin (theta (i) *pi/180) / . 
(60*cos(cglat*pi/180) ) ) ; 
end 

if  num_of_beams  ==  24 

plot (circle_lon, circle_lat , 'y* ' ) 
plot ( [ -eglon  circle_lon ( 13 ) ] , ... 
[eglat  circle_lat (13) ] , '--y' ) 
circle_lat  =  [eglat  circle_lat (1 : 24)  eglat]; 
circle_lon  =  [-eglon  circle_lon (1 : 24)  -eglon]; 
plot (circle_lon, circle_lat, 'y' ) 
else 

plot (circle_lon, circle_lat, 'y* ' ) 

plot ( [-eglon  circle_lon( (num_of_beams/2)  +  1)],... 
[eglat  circle_lat ( (num_of_beams/2)  +  l)],'--y') 
plot (circle_lon, circle_lat, 'y ' ) 
end 
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function  [cout]  =  lwad_contour ( f name, interval , overlap) 

%LWAD_CONTOUR 

%  Create  a  sediment  contour  based  on  the  selected  ping. 

% 

%  This  plot  is  created  by  determining  the  difference 

%  between  area  averaged  slope  data  and  slope  of 

%  individual  beam  data. 

% 

%  If  no  filename,  interval,  or  overlap  are  specified,  a 

%  gui  will  open  to  allow  the  user  to  select  a  file. 

%  The  interval  and  overlap  will  default  to  1  sec  and 

%  .75.   Modify  this  file  in  order  to  change  the 

%  defaults. 

% 

%  The  plot  shows  the  eg  current  position  (black  circle) , 

%  target  current  position  (red  circle) ,  and  exercise 

%  area  (black  dashed  lines) . 

% 

%  This  function  uses  code  from  the  Matlab  function 

%  ' contour. m' 

% 

%  See  also  LWAD_CONTOURF ,  CONTOUR,  LWAD. 


% 


Created:  May  15,  1999 
David  Schalm 

Last  Modified:  July  26,  1999 
David  Schalm 


if  nargin  ==  0 

load  data_drive 

[filename, pathname]  =  uigetf ile ( [data_drive  '*.*']); 

fname  =  [pathname  filename] ; 

data_drive  =  pathname; 

save  'e:\schalm\lwad  99-l\data_drive.mat '  data_drive 

%%%%%%  Default  values  %%%%% 

interval  =  1 ; 

overlap  =  .75; 

spacing  =  interval  -  interval*overlap; 

% interval  =  input ('What  is  the  interval  (.5,1,2):  '  ) ; 

%overlap  =  input ('What  is  the  overlap  (0,. 25,. 5,. 75):  '); 
5,9.9.3.9.5.3.5.5.9.9.9.0^9.5.9.9.9.5.9,9.9.5.9.9.3.0. 

else 

if  fname ( (length (fname ) -2) : length (fname ) )  ==  'mat' 

filename  =  fname ( (length (fname ) -23 ): length (fname) ) ; 
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else 

filename  =  f name ( (length (f name) -19) : length (f name) ) ; 
end 
end 

%  Load  the  data  file 

load ( f name , ' -mat ' ) 

[cglat,cglon, tgtlat, tgtlon]  =  lwad_posit (fname) ; 

[course, speed, ctrTbrng, bo tdepth, deangle]  =  lwad_prm ( fname) ; 

if  filename(5:6)  ==  'SD'  |  filename (5 : 6)  ==  'sd' 

file_type  =  filename (5 : 6) ; 
else 

file_type  =  filename (19 : 20) ; 
end 

%  Load  the  average  statistics 

area_stats_fname  =  ['e:\schalm\lwad  99-l\area_stats\ '  ... 

file_type  '_stats'  num2str (interval)  '.mat']; 
load (area_stats_f name) 
if  overlap  ==  0 

area_stats  =  stats_0; 
elseif  overlap  ==  .25 

area_stats  =  stats_25; 
elseif  overlap  ==  .5 

area_stats  =  stats_5; 
else 

area_stats  =  stats_75; 
end 

%  Create  time\ range  vector 

sound_speed  =  .81;  %nm/s 

N  =  round (interval* samplerate ) ; 

data_start  =  1; 

data_duration  =  length(bmf) /samplerate  -  (1/ samplerate) . . . 

+  data_start; 
time  =  (data_start : 1/samplerate :data_duration) ' ; 
[mtime,ao]  =  mmean (time, N, overlap) ; 
mrange  =  sound_speed*mtime/2 ; 

%  Create  bearing  vector 

num_of_beams  =  size (bmf , 2 ) ; 

theta_one  =  ctrTbrng  -  num_of_beams*5/2  +  5/2; 

theta  =  theta_one: 5 : ( theta_one+ (num_of_beams-l) *5) ; 

index  =  0 ; 

for  i  =  l:l:num  of  beams 
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if  sum(bmf ( : , i)  ==0)  >  0 

%  skip  this  beam 
else 

log_data  =  2  0*logl0 (abs (bmf ( : , i) ) ) ; 

index  =  index  +  1; 

thetal (index)  =  theta(i); 

mean_data  =  mmean (log_dat a, N, overlap) ; 

slope_data  =  (gradient (mean_data) ) ; 

diff_data (index, : )  =  slope_data  -  area_stats (1, : ) ; 
end 
end 

cl  =  contourc (mrange, thetal, diff_data) ; 

i  =  1; 

while  i  <  length (cl) 
nexti  =  i  +  cl(2,i); 

c_new ( : , i )  =  c  1  (  :  ,  i )  ; 

%  convert  range /brng  to  lat/lon 
for  k  =  (i+1) : 1 : nexti 

c_new(l,k)  =  -(cglon  -  sin (cl (2 , k) *pi/180) * . . . 

cl (l,k) / (60*cos (cglat*pi/180) ) ) ; 
c_new(2,k)  =  cglat  +  cos (cl (2,k) *pi/180) * . . . 
cl(l,k) /60; 
end 

i  =  nexti  +  1; 
end 
c  =  c_new ; 

figure 

lims  =  [-84.4  -84  25.5  25.8]; 

colortab  =  get (newplot, ' colororder' ) ; 
[mc,nc]  =  size (colortab) ; 

%  Create  the  plot 

newplot; 

if  -ishold 

view (3);  grid  on 

set (gca, 'xlim' , lims (1:2) , 'ylim' , lims (3:4) ) 
end 

limit  =  size (c, 2); 
i  =  1; 
h  =  [  ]  ; 
color_h  =  [ ] ; 
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while (i  <  limit) 

z_level  =  c (1, i)  ; 
npoints  =  c(2,i); 
nexti  =  i+npoints+1; 

xdata  =  c (1, i+1 : i+npoints) ; 
ydata  =  c (2 , i+1 : i+npoints) ; 

zdata  =  z_level  +  0*xdata;   %  Make  zdata  the  same  size 

%  as  xdata 

%  Create  the  patches  or  lines 

cu  =  patch ( 'XData' , [xdata  NaN] , 'YData' , [ydata  NaN] ,  . . 

'ZData' , [zdata  NaN] , 'CData' , [zdata  NaN] ,  ... 

' f acecolor ' , ' none ' , ' edgecolor ' , ' flat ' , . . . 

'userdata' , z_level) ; 
h  =  [h;  cu(:)] ; 

color_h  =  [color_h  ;  z_level] ; 
i  =  nexti; 
end 

for  i  =  1: length (h) 

set (h(i) , 'Zdata' ,  []  )  ; 
end 

view (2 ) ; 

set (gca, 'box' , 'on') ; 

axis  equal,  grid  on,  zoom  on,  hold  on 

axis([-84.5  -83.9  25.4  25.9]) 

xlabelC Longitude  (W)  '  )  ,  ylabel  ( 'Latitude  (N)  '  ) 

warning  off 

title( ['Event  31/'  filename (13 : 16)  ',  '  ... 

ping_date  ' ,  '  ping_time] ) 
caxis([-12  12]) 
colorbar 

plot (-cglon, cglat, 'ko' ) 
plot (-tgtlon, tgtlat , ' ro' ) 

%  Plot  eg  box 

plot([-84.27  -84.27  -84.11  -84.11  -84.27],... 

[25.61  25.71  25.71  25.61  25.61], 'k— -') 
zoom  on 
brighten (1) 

if  nargout  >  0 

cout  =  c ; 
end 
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function  [cout]  =  lwad_contourf ( f name, interval, overlap) 

% LWAD_CONTOURF 

%  Create  a  filled  sediment  contour  based  on  the  selected 

%  ping. 

% 

%  This  plot  is  created  by  determining  the  difference 

%  between  area  averaged  slope  data  and  slope  of 

%  individual  beam  data. 

% 

%  If  no  filename,  interval,  or  overlap  are  specified,  a 

%  gui  will  open  to  allow  the  user  to  select  a  file. 

%  The  interval  and  overlap  will  default  to  1  sec  and 

%  .75.   Modify  this  file  in  order  to  change  the 

%  defaults. 

% 

%  The  plot  shows  the  eg  current  position  (black  circle) , 

%  target  current  position  (red  circle) ,  and  exercise 

%  area  (black  dashed  lines) . 

% 

%  This  function  uses  code  from  the  Matlab  function 

%  'contourf.m' 

% 

%  See  also  LWAD_CONTOUR ,  CONTOURF,  LWAD. 

%  Created:  May  15,  1999 

%         David  Schalm 

% 

%  Last  Modified:  July  26,  1999 

%         David  Schalm 

if  nargin  ==  0 

load  data_drive 

[filename, pathname]  =  uigetf ile ( [data_drive  '*.*']); 

fname  =  [pathname  filename] ; 

data_drive  =  pathname; 

save  'e:\schalm\lwad  99-1 \data_drive. mat '  data_drive 

%%%%%%  Default  values  %%%%% 

interval  =  1 ; 

overlap  =  .75; 

% interval  =  input ('What  is  the  interval  (.5,1,2):  '); 

%overlap  =  input ('What  is  the  overlap  (0,. 25,. 5,. 75):  '); 


else 

if  fname ( (length (fname) -2) : length (fname) )  ==  'mat' 

filename  =  fname ( (length (fname ) -23) : length (fname) ) ; 
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else 

filename  =  fname ( (length ( fname) -19) : length (fname) ) ; 
end 
end 

%  Load  the  data  file 

load ( f name , ' -mat ' ) 

[cglat, cglon, tgtlat, tgtlon]  =  lwad_posit (fname) ; 

[course, speed, ctrTbrng, bo tdepth, deangle]  =  lwad_prm ( fname ] 

spacing  =  interval  -  interval * over lap; 


if  filename (5: 6)  ==  'SD'    filename (5 : 6)  ==  'sd' 

file_type  =  filename (5 : 6) ; 
else 

file_type  =  filename (19 : 20) ; 
end 

%  Load  the  average  statistics 

area_stats_fname  =  ['e:\schalm\lwad  99-l\area_stats\ '  . 

file__type  '_stats'  num2str (interval)  '.mat']; 
load (area_s tat s_f name) 
if  overlap  ==  0 

area_stats  =  stats_0; 
elseif  overlap  ==  .25 

area_stats  =  stats_25; 
elseif  overlap  ==  .5 

area_stats  =  stats_5; 
else 

area_stats  =  stats_75; 
end 

%  Create  time\range  vector 

sound_speed  =  .81;  %nm/s 

N  =  round (interval* samplerate ) ; 

data_start  =  1; 

data_duration  =  length (bmf) /samplerate  -  (1 /samplerate) 

+  data_start; 
time  =  (data_start : 1 / samplerate :data_durat ion) ' ; 
[mtime,ao]  =  mmean (time, N, overlap) ; 
mrange  =  sound_speed*mtime/2 ; 

%  Create  bearing  vector 

num_of_beams  =  size (bmf, 2); 

theta_one  =  ctrTbrng  -  num_o f _beams * 5 / 2  +  5/2; 

theta  =  theta_one: 5 : (theta_one+ (num_of_beams-l) *5) ; 
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index  =  0 ; 

for  i  =  1 : 1  :num_of__beams 

if  sum(bmf(:,i)  ==0)  >  0 

%  skip  beams  w/  bad  data 
else 

log_data  =  20*logl0 (abs (bmf ( : , i) ) ) ; 
index  =  index  +  1 ; 
thetal (index)  =  theta(i); 
mean_data  =  mmean (log_data, N, overlap) ; 
slope_data  =  (gradient (mean_data) ) ; 

diff_data (index, : )  =  (slope_data  -  area_stats (1, : ) ) ; 
end 
end 

diff_data  =  diff_data/ spacing; 


/  / 
/  / 


lin  = 

col  = 

x  =  mrange; 

y  =  thetal ; 

z  =  diff_data; 

nv  =  [ ] ; 

%nv  =  4; 

if  (size(y, 1) ==1) ,  y=y' ;  end; 

if  (size (x, 2) ==1) ,  x=x' ;  end; 

[mz, nz]  =  size (z)  ; 

lims  =  [-84.4  -84  25.5  25.8]; 


i  =  f ind(isf inite (z) ) ; 
minz  =  min (z (i) )  ; 
maxz  =  max ( z ( i ) ) ; 


if  length (nv)  <=  1 
if  isempty(nv) 

CS=contourc ( [minz  maxz 
else 

CS=contourc ( [minz  maxz 
end 


minz  maxz] ) ; 
minz  maxz] ,nv) ; 


%  Find  the  levels 

ii  =  1; 

nv  =  minz;  %  Include  minz  so  that  the  contours  are 

%totally  filled 

while  (ii  <  size(CS,2)), 

nv= [nv  CS(1, ii) ] ; 

ii  =  ii  +  CS(2,ii)  +  1; 
end 
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end 

%  Don't  fill  contours  below  the  lowest  level  specified  in 

%  nv.   To  fill  all  contours,  specify  a  value  of  nv  lower 

%  than  the  minimum  of  the  surface. 

draw_min=0 ; 

if  any(nv  <=  minz) , 

draw_min=l; 
end 

%  Get  the  unique  levels 
nv  =  sort ( [minz  nv (:)']); 
zi  =  [1,  find(diff (nv) )+l]  ; 
nv  =  n v ( z  i ) ; 

%  Surround  the  matrix  by  a  very  low  region  to  get  closed 
%  contours,  and  replace  any  NaN  with  low  numbers  as  well. 

zz= [  repmat (NaN, l,nz+2)  ;  repmat (NaN,mz, 1)  z  ... 

repmat (NaN,mz, 1)  ;  repmat (NaN, 1, nz+2 )] ; 
kk=f ind(isnan (zz ( : ) ) ) ; 
zz (kk) =minz-le4* (maxz-minz) +zeros (size (kk) ) ; 

xx  =  [2*x ( : , 1) -x( : , 2 ) ,  x,  2*x ( : , nz) -x ( : , nz-1 ) ] ; 
yy  =  [2*y(l, :)-y(2, :) ;  y;  2*y (mz, : ) -y (mz-1, : ) ] ; 
if  (min(size (yy) ) ==1) , 

[CS,msg] =contours (xx,yy, zz,nv) ; 
else 

[CS,msg] =contours (xx( [  1  l:mz  mz] , : ) , yy ( : , [1  l:nz  nz] ) . 
, zz,nv) ; 
end; 
if  -isempty  (msg) ,  error(msg);  end 

c_new  =  [ ] ; 

i  =  1; 

while  i  <  length (CS) 

nexti  =  i  +  CS(2,i) ; 

c_new(:,i)  =  CS(:,i); 

%   convert   range/brng   to    lat/lon 
for  k  =    (i+1)  :  lmexti 

c_new(l,k)    =    - (cglon   -    sin(CS(2,k) *pi/180) * . . . 

CS(l,k) / (60*cos (cglat*pi/180) ) ) ; 
c_new(2,k)  =  cglat  +  cos (CS (2 ,k) *pi/180) * . . . 
CS(l,k) /60; 
end 
i  =  nexti  +  1; 
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end 

CS  =  c_new; 

figure 

%  Find  the  indices  of  the  curves  in  the  c  matrix,  and  get 

%  the  area  of  closed  curves  in  order  to  draw  patches 

%  correctly. 

ii  =  1; 

ncurves  =  0 ; 

I  =  []; 

Area= [ ] ; 

while  (ii  <  size(CS,2)), 

nl=CS(2, ii) ; 

ncurves  =  ncurves  +  1; 

I (ncurves)  =  ii; 

xp=CS(l, ii+ (l:nl) ) ;   %  First  patch 

yp=CS(2,ii+(l:nl) ) ; 

Area (ncurves) = sum (  dif f (xp) . * (yp (1 :nl-l) +yp (2 :nl) ) /2  ) ; 

ii  =  ii  +  nl  +  1; 
end 

newplot ; 
if  -ishold, 

view (2 ) ; 

set (gca, 'xlim' , lims (1:2) , 'ylim' , lims (3:4) ) 
end 

%  Plot  patches  in  order  of  decreasing  size.  This  makes  sure 
%  that  all  the  leves  get  drawn,  not  matter  if  we  are  going 
%  up  a  hill  or  down  into  a  hole.  When  going  down  we  shift 
%  levels  though,  you  can  tell  whether  we  are  going  up  or 
%  down  by  checking  the  sign  of  the  area  (since  curves  are 
%  oriented  so  that  the  high  side  is  always  the  same  side) . 
%  Lowest  curve  is  largest  and  encloses  higher  data  always. 

H=  [  ]  ; 

[FA, IA] =sort (-abs (Area) ) ; 
if  -isstr (get (gca, ' color ')) , 

bg  =  get (gca, 'color ') ; 
else 

bg  =  get (gcf, ' color ') ; 
end 
if  isempty(col) 

edgec  =  get (gcf , 'defaultsurfaceedgecolor ' ) ; 
else 

edgec  =  col ; 
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end 

if  isempty (lin) 

edgestyle  =  get (gcf , 'defaultpatchlinestyle' ) ; 
else 

edgestyle  =  lin; 
end 
for  jj=IA, 

nl=CS(2,I(jj))  ; 
lev=CS(l,I(jj) )  ; 
if  (lev  ~=minz   draw_min  )  , 
xp=CS(l,I(jj)  +  (l:nl))  ; 
yp=CS(2/I(jj)  +  (l:nl) )  ; 

if  (sign (Area (jj ) )  ~=sign (Area (IA(1) ) )  ), 
kk=f ind(nv==lev) ; 
if  (kk>l+sum(nv<=minz) * (~draw_min) ) , 

lev=nv(kk-l) ; 
else 

lev=NaN;  %  missing  data  section 

end; 
end; 

if  (isf inite (lev) ) , 

H= [H;patch (xp,yp, lev, ' facecolor ' , ' flat' , . . . 

'edgecolor' ,edgec, ' lines tyle' , edgestyle) ] ; 
else 

H= [H;patch(xp,yp, lev, ' facecolor' ,bg, ... 

' edgecolor' , edgec, ' lines tyle' , edgestyle) ] ; 
end; 
end; 
end; 

numPatches  =  length (H) ; 
if  numPatches>l 

for  i=l : numPatches 

set(H(i),  ' faceof f setfactor ' ,  0,  ... 

' faceof fsetbias ' ,  (le-3) + ( numPatches- i) / . . . 
( numPatches- 1) /3  0)  ; 
end 
end 

%  Set  the  plot  options 

axis  equal,  grid  on,  zoom  on,  hold  on 

axis([-84.5  -83.9  25.4  25.9]) 

xlabelC Longitude  (W) ' ) ,  ylabel (' Latitude  (N) ' ) 

warning  off 

title ( [ 'Ping: '  filename (13 : 16)  ... 

',Mode:  SD,  Interval:  1,  Overlap:  .75']) 
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%  Uncomment  the  next  line  for  a  b/w  plot 
%colormap (gray) 

caxis( [-12  12]  ) 

colorbar 

plot (-cglon, cglat , 'ko' ) 

plot (-tgtlon, tgtlat, 'ro' ) 

%  Plot  eg  box 

plot([-84.27  -84.27  -84.11  -84.11  -84.27], 

[25.61  25.71  25.71  25.61  25 . 61] , ' k-- ' ) 
zoom  on 

if  nargout  >  0 

cout  =  c; 
end 
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%  COMPOS IT 

% 

%  This  script  file  will  create  a  composite  graph  of  many 

%  separate  pings.   For  now,  it  will  only  work  for  SD  mode 

%  pings . 

%  Location  of  the  pings  to  combine 

filepath  =  {'e:\schalm\LWAD  99-1  Data\Tape  17\sd\'}; 

%%%%  Default  values  (Modify  if  desired)  %%%% 

interval  =  1; 

overlap  =  .75; 

spacing  =  interval  -  interval*overlap; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Load  the  area  stats  file 

area_stats_fname  =  ['e:\schalm\lwad  99-l\area_stats\ ' . . . 
'sd_stats'  num2str (interval)  '.mat']; 

load (area_stats_f name) 

if  overlap  ==  0 

area_stats  =  stats_0; 
elseif  overlap  ==  .25 

area_stats  =  stats_25; 
elseif  overlap  ==  . 5 

area_stats  =  stats_5; 
elseif  overlap  ==  .75 

area_stats  =  stats_75; 
end 

%  Get  a  directory  of  the  pings  to  combine 

D  =  dir (char (filepath) ) ; 

[fname{l : length (D) } ]  =  deal (D. name) ; 

%  get  rid  of  ' . '  and  ' . . '  filenames 
f name  =  f name ( 3 : length ( f name ) ) ; 

%  Select  the  Lat/Lon  you  want  to  look  at 
%  and  the  grid  resolution 
Ion  =  -84.45: .01:-83.9; 
lat  =  25.45: .01:25.9; 

compos it_graph  =  zeros (length( Ion) , length( lat) )  ; 
ref_matrix  =  zeros (length (Ion) , length (lat) ) ; 

%  Fill  in  the  values  of  compos it_graph  one  ping  at  a  time 
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for  m  =  1:1:1 

%for  m  =  1 : 1 : length (f name) 

load( [char (f ilepath)  char (fname (m) ) ] , ' -mat ' ) 

disp ([ 'Current  file  is  '  char (f ilepath)  char (fname {m} )] ) 

[cglat , cglon, tgtlat, tgtlon]  =  lwad_posit . . . 

( [char ( f ilepath)  char ( fname (m) )  ]  )  ; 
[course, speed, ctrTbrng, botdepth, deangle]  =  ... 

lwad_prm( [char (f ilepath)  char (fname (m) )  ]  )  ; 
cglon  =  -cglon; 

for  lon_index  =  1 : 1 : length (Ion) 

for  lat_index  =  1 : 1 : length (lat) 

%  check  to  see  if  lat/lon  is  in  range  of  the  eg 
dx  =  60*cos (cglat*pi/180) * (lon(lon_index) -cglon) ; 
dy  =  60* (lat (lat_index)  -  cglat); 
range  =  sqrt(dx^2  +  dy/"2); 

if  range  >  9.9 

%  do  nothing  because  this  is  out  of  range 

%  of  the  eg 
else 

brng  =  atan2 (dy, dx) ; 

%  Convert  bearings  to  true  (north  up)  bearings 
if  brng  <  0 

brngl  =  -1* (brng*180/pi-90) ; 
else 

brngl  =  90-brng*180/pi; 
end 

if  brngl  <  0 

brngl  =  brngl+3  60; 
end 
brngl; 
beam_num  =  ceil (brngl/5) ; 

%  Skip  over  data  with  0  values  -->  cause  errors 
%  in  statistics 

if  sum(bmf ( : , beam_num)  ==0)  >  0 
%do  nothing 

%  Skip  over  beams  that  are  +/-  45  degrees 

%  of  the  stern 
elseif  course  >  075  &  course  <  105  &  ... 
beam_num  >  45  &  beam_num  <  63 

%  do  nothing 
elseif  course  >  165  &  course  <  195  &  ... 
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( beam_num  <  9   beam_num  >  63) 
%  do  nothing 
else 

N  =  round ( interval *samplerate) ; 
log_data  =  20*logl0 (abs (bmf ( : ,beam_num) ) ) ; 
mean_data  =  mmean (log_data,N, overlap) ; 
slope_data  =  (gradient (mean_data) ) ; 
diff_data  =  (slope_data  -  area_stats . . . 

(1, : ) ) /spacing; 
sound_speed  =  .81;  %nm/s 
data_time  =  range*2/sound_speed; 
index  =  round ( (data_time- .25) / . 25) ; 
if  index  <  1 
index  =  1 ; 
end 

diff_point  =  dif f_data (index) ; 
composit_graph(lon_index, lat_index)  =  ... 
composit_graph(lon_index, lat_index) + . . . 
dif f_point ; 
ref_matrix(lon_index, lat_index)  =  ... 
ref_matrix (lon_index, lat_index)  +  1; 
end 
end 
end 
end 
end 

%  Divide  the  lat/lon  matrix  values  by  the  number  of 

%  entries  at  each  point 

composit_graph  =  composit_graph. /ref_matrix; 

%  Plot  the  composite  graph 

figure 

contourf (Ion, lat, compos it_graph' ) ,  hold  on 

%  Uncomment  the  next  line  for  a  b/w  plot 
%colormap (gray) 
caxis( [-8  8] ) , 
grid  on,colorbar 

xlabel ( ' Longitude  (W) ' ) ,  ylabel ( ' Latitude  (N) ' ) 
title ([ 'Composite  Plot']) 

%  Plot  eg  box 

plot([-84.27  -84.27  -84.11  -84.11  -84.27],... 

[25.61  25.71  25.71  25.61  25.61], 'k— ') 
zoom  on 
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%LWAD_FILE 

% 

%   Use  this  script  to  open  and  load  an  LWAD  99-1  file 

load  data_drive 

[filename, pathname]  =  uigetf ile ( [data_drive] ) ; 

fname  =  [pathname  filename] ; 

data_drive  =  pathname; 

save  'e:\schalm\lwad  99-l\data_drive.mat '  data_drive 

%  Load  the  data  file 
load ( fname ,  ' -mat ' ) 
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